薬学のための医療統計学

#薬学のための医療統計学　

library(dplyr)

#検定


#2群間の平均値の差の検定

#対応がある場合<=>1対の標本から得られる

#対応がない場合<=>独立な標本から得られる

#正規分布している<=>パラメトリックな検定

#正規分布していない<=>ノンパラメトリックな検定

#(1)対応のある２群間のt-test(paired t-test:パラメトリックな検定):

#H：各個体において投与前と投与後の値は等しい

x=c(38.6,38.5,39,38.3,38.7,38.1,38.5,38.2,38.4,38.8)

y=c(37.2,37.6,38.5,37.5,38.1,36.8,37.8,37.2,37.6,38.2)

Z=x-y;Z_bar=mean(Z);SD=sd(Z);n=length(x)

t_cal=abs(Z_bar)/(SD/sqrt(n))

qt(1-0.05/2,df=n-1)

#(2)対応のある場合の2群間のノンパラメトリック検定法（Wilcoxon検定):

#投与前
x=c(8,6,3.5,2,7,6.5,5,8.5,5.5,7.5)

#投与後
y=c(7.5,4.3,3.5,2.4,6.4,3.2,4.8,8,5.8,5.3)

D=-(x-y);

D_dat=data.frame(D=(D))%>%mutate(ab_D=abs(D))%>%filter(ab_D>0)

D_dat=D_dat[order(D_dat$ab_D),]

D_dat=data.frame(no=c(1:nrow(D_dat)),D_dat)

or_D_dat=D_dat%>%group_by(ab_D)%>%summarise(D=mean(D),no=mean(no),n=n())%>%mutate(no2=no^2)

N=sum(or_D_dat$n)

option=1

if(option==1){

#同順位の場合

R0=N*(N+1)/4-qnorm(1-0.05/2)*sqrt(sum(or_D_dat$no2)/4)-0.5

}else{

#大標本の場合

R0=N*(N+1)/4-qnorm(1-0.05/2)/sqrt(N*(N+1)*(2*N+1)/24)-0.5

}

R_plus=sum(or_D_dat$no[or_D_dat$D>0]);R_minus=sum(or_D_dat$no[or_D_dat$D<0])

R_cal=min(R_plus,R_minus)

#R_cal<R0なら棄却

R_cal<R0


#対応のない場合

#等分散性の検定（例題６．５）：

x=c(122,81,93,150,145,84,106,125,136,115)

y=c(130,125,110,173,152,186,108,165,157)

Sx=sum((x-mean(x))^2)/(length(x)-1)

Sy=sum((y-mean(y))^2)/(length(y)-1)

f=Sy/Sx

qf(1-0.05/2,length(y)-1,length(x)-1)

#有意差あり＝＞Welch-test;有意差なし＝＞Student-t-test



#student-t-test(パラメトリックかつ等分散）：

x=c(122,81,93,150,145,84,106,125,136,115)

y=c(130,125,110,173,152,186,108,165,157)

ave_x=mean(x);ave_y=mean(y)

Sx=sum((x-mean(x))^2)

Sy=sum((y-mean(y))^2)

t_cal=abs(ave_x-ave_y)/sqrt((1/length(x)+1/length(y))*((Sx)+(Sy))/(length(x)+length(y)-2))

qt(1-0.05/2,length(x)+length(y)-2)



#Welch-test(パラメトリックかつ非等分散性):

x=c(122,81,93,150,145,84,106,125,136,115)

y=c(140,151,145,173,152,164,168,171,157)

ave_x=mean(x);ave_y=mean(y)

Sx2=sum((x-mean(x))^2)/(length(x)-1)

Sy2=sum((y-mean(y))^2)/(length(y)-1)

t_cal=abs(ave_x-ave_y)/sqrt(Sx2/length(x)+Sy2/length(y))

pthi=(Sx2/length(x)+Sy2/length(y))^2/((Sx2/length(x))^2/(length(x)-1)+(Sy2/length(y))^2/(length(y)-1))

if(pthi==floor(pthi)){

bound=qt(1-0.05/2,df=pthi)

}else{

bound=qt(1-0.05/2,df=length(y))*(120/length(x)-120/pthi)/(120/length(x)-120/length(y))+qt(1-0.05/2,df=length(x))*(120/pthi-120/length(y))/(120/length(x)-120/length(y))

}


#Mann-Whitney U-test(ノンパラメトリック）：

A=c(2,5,2,3,4,3,3);B=c(7,6,5,5,8,5)

vec=c(A,B)

vec_dat=data.frame(no=c(1:length(vec)),vec=sort(vec))

vec_dat2=vec_dat%>%group_by(vec)%>%summarise(no=mean(no),n=n())

or_A=or_vec[1:length(A)];or_B=or_vec[-c(1:length(A))]

vec_order=c()

for(j in 1:length(vec)){

vec_order=c(vec_order,vec_dat2$no[vec_dat2$vec==vec[j]])

}

R_A=sum(vec_order[c(1:length(A))])

R_B=sum(vec_order[-c(1:length(A))])

U_A=length(A)*length(B)+length(B)*(length(B)+1)/2-R_B

U_B=length(A)*length(B)+length(A)*(length(A)+1)/2-R_A

U_cal=min(U_A,U_B)


#多重比較法

#tukey-kramer 

data=data.frame(experience_times=c(1,2,3),group1=c(13,17,13),group2=c(20,19,17),group3=c(12,15,15),group4=c(22,19,18))

x_ave=apply(data[,-1],2,mean)

V=apply(data[,-1],2,var)

pthi_E=prod(dim(data[,-1]))-ncol(data[,-1])

V_E=sum((nrow(data)-1)*V)/pthi_E

t=array(0,dim=c(ncol(data[,-1]),ncol(data[,-1])))

for(i in 1:ncol(t)){
for(j in 1:nrow(t)){

t[i,j]=(x_ave[i]-x_ave[j])/sqrt(V_E*(1/nrow(data)+1/nrow(data)))

}}

#Dunnet

data=data.frame(patients=c(1,2,3),group1=c(110,100,105),group2=c(95,100,95),group3=c(120,100,105),group4=c(91,90,89))

x_ave=apply(data[,-1],2,mean)

V=apply(data[,-1],2,var)

pthi_E=prod(dim(data[,-1]))-ncol(data[,-1])

V_E=sum((nrow(data)-1)*V)/pthi_E

t=array(0,dim=c(ncol(data[,-1]),ncol(data[,-1])))

for(j in 1:nrow(t)){

t[1,j]=(x_ave[1]-x_ave[j])/sqrt(V_E*(1/nrow(data)+1/nrow(data)))

}