医学・薬学データの統計解析　データの整理から交互作用多重比較まで　廣津千尋先生

```{r}
#薬の検定

n0=100;n1=200

#既存薬のデータ

y0=rnorm(n0,log(5),1)

#非検薬のデータ

y1=rnorm(n1,log(6),1)

x0=(y0);x1=(y1)

sigma_hat=(sum((x0-mean(x0))^2)+sum((x1-mean(x1))^2))/(n0+n1-2)

print(paste0("対数を取った場合の90%信頼区間は",mean(x0)-mean(x1)+sqrt(1/n0+1/n1)*sigma_hat*qt(0.05,n0+n1-2),",",mean(x0)-mean(x1)+sqrt(1/n0+1/n1)*sigma_hat*qt(1-0.05,n0+n1-2)))

#90%信頼区間が0.8～1.25に収まっていることが望ましい

print(paste0("90%信頼区間は",exp(mean(x0)-mean(x1)+sqrt(1/n0+1/n1)*sigma_hat*qt(0.05,n0+n1-2)),",",exp(mean(x0)-mean(x1)+sqrt(1/n0+1/n1)*sigma_hat*qt(1-0.05,n0+n1-2))))

#推定値は５％信頼区間に含まれるため、棄却されない（有意差がない）

mean(x1)-mean(x0)





```

```{r}

#2薬剤の優劣性の検定(Wilcoxson)　p.70~

data=data.frame(medicine=c("AMPC","S6472"),無効=c(3,8),やや無効=c(8,9),有効=c(30,29),著効=c(22,11),合計=c(63,57))

C=apply(data[,2:(ncol(data)-1)],2,sum)

w=c()

for(j in 1:4){
  
c_sub=C

c_sub[j]=(c_sub[j]+1)/2
  
w=c(w,sum(c_sub[1:j]))  
  
}

n1=data$合計[1];n2=data$合計[2];n=n1+n2

z1=c(data[1,2:5])

W1=sum(w*as.numeric(z1))

E_W1=(n1/n)*sum(w*C)

V_W1=(n1*n2)*(sum(C*w^2)-sum(C*w)^2/n)/(n*(n-1))

sigma_w=sum(C*w^2)/n-(sum(C*w)/n)^2

W=sqrt((n-1)/(n1*n2*sigma_w))*(W1-E_W1)

qnorm(1-0.05/2)

1-pnorm(W)

```

```{r}

data=data.frame(mg25=c(1.55,1.63,1.49,1.53,2.14),mg200=c(1.78,1.93,1.8,2.07,1.70))

n=prod(dim(data))

E_X=sum(as.matrix(data))/n

V_X=sum(as.matrix(data)^2)/n-E_X^2

n1=length(data$mg25);n2=length(data$mg200)

COV_X=-n1*V_X/(n2*(n-1))

E_X_bar=E_X

V_X_bar=V_X/(n-1)

y1=data$mg25;y2=data$mg200

u=(mean(y2)-E_X)/sqrt(V_X_bar)

#標準正規分布

qnorm(1-0.05/2)

```

```{r}

#累積χ二乗法

dt=data.frame(y=c(y1,y2),order=order(c(y1,y2)),rank=c(rep(1,5),rep(2,5)))

dt=dt[order(dt$order),]

dt=data.frame(mg200=ifelse(dt$rank==2,1,0),mg25=ifelse(dt$rank==1,1,0))

n=apply(dt,2,sum)

J=4;

mu=J/sum(n)

mat=cbind(t(t(apply(dt[1:J,],2,sum))),t(t(apply(dt[(J+1):nrow(dt),],2,sum))))

sig2=J*(sum(n)-J)/(sum(n)^2)

E_W1=n[1]*(sum(n)-J)/sum(n)

V_W1=(n[1]*n[2]*J*(sum(n)-J))/((sum(n)-1)*sum(n)^2)

W1=-(mat[1,1]-n[1]*J/sum(n))/sqrt((n[1]*n[2]*J*(sum(n)-J))/((sum(n)-1)*sum(n)^2))

kai=W1^2

((sum(n)-1)*(sum(n)*mat[1,1]-n[1]*J)^2)/(n[1]*n[2]*J*(sum(n)-J))

#χ二乗分布(df=1)

qchisq(1-0.05/2,df=1)

```

```{r}

#ハンディキャップ検定

n1=100;n2=200;

y1=rpois(1,lambda=sample(c(1:100),n1))

y2=rpois(1,lambda=sample(c(1:100),n2,replace=T))

abs(y1/n1-y2/n2)/sqrt(((y1+y2)/(n1+n2))*(1-(y1+y2)/(n1+n2))*(1/n1+1/n2))

qnorm(1-0.05/2)

#非劣性検証を仮定(p1>p2)

delta=0.5;p2=0.3

p1=p2+delta;

p1-p2

V1=p1*(1-p1)/n1+(p1-delta)*(1-p1+delta)/n2

U=(y1/n1-y2/n2-(p1-p2))/sqrt(V1)

#棄却(TRUE)

abs(U)>qnorm(1-0.05/2)

alpha=0.05;beta=0.95

V0=(1/n1+1/n2)*((y1+y2)/(n1+n2))*(1-(y1+y2)/(n1+n2))

p=sqrt(V0/V1)

#必要サンプルサイズ(n1=n2)

(p1*(1-p1)+(p1-delta)*(1-p1+delta))*((p*qnorm(1-alpha/2)+qnorm(beta))/(p1-p2))^2

print(paste0("検出力：",1-beta))

```

```{r}

#ハンディキャップ検定２

n1=100;n2=200;

y1=rpois(n1,lambda=sample(c(1:100),n1))

y2=rpois(n2,lambda=sample(c(1:100),n2,replace=T))

delta=2;alpha=0.05;beta=0.05

t=(mean(y1)+delta-mean(y2))/sqrt((1/n1+1/n2)*(sum((y1-mean(y1))^2)+sum((y2-mean(y2))^2))/(n1+n2-2))

qt(1-alpha/2,df=n1+n2-2)

#サンプルサイズ(検出力が1-betaとなる場合)

mu1=mean(y1);mu2=mean(y2)

fun=function(m){return(2*(qt(1-alpha,df=2*(m-1))+qt(1-beta,df=2*(m-1)))^2/((delta+mu1-mu2)/sqrt((sum((y1-mean(y1))^2)+sum((y2-mean(y2))^2))/(n1+n2-2)))^2-m)}

times=100
newton_raphson=data.frame(times=1:times,value=0,f=0,df=0)
value=20;h=0.001
for(j in 1:nrow(newton_raphson)){
f<-fun(value)
df<-(fun(value+h)-fun(value))/h

  newton_raphson$f[j]=f
  newton_raphson$df[j]=df
  value=-f/df+value

  newton_raphson$value[j]=value
}

solution=newton_raphson$value[nrow(newton_raphson)]

print(paste0("必要サンプルサイズ：",solution))

```
```{r}

#Dunnet p.91

data=data.frame(medicines=c("AMPC","S6472","CCL"),無効=c(3,8,2),やや有効=c(8,9,11),有効=c(30,29,33),著効=c(22,11,17))

C=apply(data[,2:ncol(data)],2,sum)

n=apply(data[,2:ncol(data)],1,sum)

w=c()

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

c_sub=C

c_sub[j]=(c_sub[j]+1)/2

w=c(w,sum(c_sub[1:j]))  

}

sigma_w=sum(C*w^2)/sum(n)-(sum(w*C)/sum(n))^2

sigma_w=sqrt(sigma_w)

wz=c()

for(j in 1:nrow(data)){
  
wz=c(wz,sum(w*data[j,2:ncol(data)]))  
  
}

W_mat=array(0,dim=c(nrow(data),nrow(data)))

for(j in 1:nrow(W_mat)){
for(i in 1:ncol(W_mat)){  
  
if(i!=j){  
  
W_mat[i,j]=sqrt((sum(n_vec)-1)/sum(n_vec))*(1/sigma_w)*(wz[i]/n[i]-wz[j]/n[j])*(1/n[i]+1/n[j])^(-1/2) 
   
}  
  
}
}

n_vec=n

```

```{r}

#Dunnet p.94

l=1

w0=-1

r=n_vec[-l]/(n_vec[-l]+n_vec[l])

func=function(w1){

return(prod(pnorm((w0-sqrt(r)*w1)/sqrt(1-r)))*dnorm(w1))
  
}

a=-100;H=abs(a)*2;n=100000

Y=(H/3)*(1/(n-1))*(func(a)+func(a+H));

for(j in 1:(n-2)){
  coefficient=ifelse(j-floor(j/2)*2==0,4,2)
  Y=Y+((H*coefficient)/(3*(n-1)))*func(a+(j*H)/(n-1))

}
#simpson積分値
print(Y)




```

```{r}

#Dunnet w0 p.94

l=1

interval=seq(-5,5,0.1)

result=data.frame(interval,f=0)

for(b in 1:length(interval)){

w0=interval[b]

r=n_vec[-l]/(n_vec[-l]+n_vec[l])

func=function(w1){

return(prod(pnorm((w0-sqrt(r)*w1)/sqrt(1-r)))*dnorm(w1))
  
}

a=-100;H=abs(a)*2;n=100000

Y=(H/3)*(1/(n-1))*(func(a)+func(a+H));

for(j in 1:(n-2)){
  coefficient=ifelse(j-floor(j/2)*2==0,4,2)
  Y=Y+((H*coefficient)/(3*(n-1)))*func(a+(j*H)/(n-1))

}
#simpson積分値
print(Y)

result$f[b]=Y

}

plot(result$interval,result$f,type="o",col=2)

```

```{r}

#tukey p.95

l=1

interval=seq(0,5,0.1)

result=data.frame(interval,f=0)

for(b in 1:length(interval)){

w0=interval[b]

r=n_vec[-l]/(n_vec[-l]+n_vec[l])

func=function(w1){

z=0  
  
z=z+prod(rep((pnorm(w1)-pnorm(w1-sqrt(2)*w0)),length(C)-1))*dnorm(w1)*length(C)

return(z)
  
}

a=-100;H=abs(a)*2;n=100000

Y=(H/3)*(1/(n-1))*(func(a)+func(a+H));

for(j in 1:(n-2)){
  coefficient=ifelse(j-floor(j/2)*2==0,4,2)
  Y=Y+((H*coefficient)/(3*(n-1)))*func(a+(j*H)/(n-1))

}
#simpson積分値
print(Y)

result$f[b]=Y

}

plot(result$interval,result$f,type="o",col=2)



```

```{r}

n0=100;n1=200

#既存薬のデータ

y0=rnorm(n0,log(5),1)

#非検薬のデータ

y1=rnorm(n1,log(6),1)

x0=(y0);x1=(y1)

sigma_hat=(sum((x0-mean(x0))^2)+sum((x1-mean(x1))^2))/(n0+n1-2)

print(paste0("対数を取った場合の90%信頼区間は",mean(x0)-mean(x1)+sqrt(1/n0+1/n1)*sigma_hat*qt(0.05,n0+n1-2),",",mean(x0)-mean(x1)+sqrt(1/n0+1/n1)*sigma_hat*qt(1-0.05,n0+n1-2)))

#90%信頼区間が0.8～1.25に収まっていることが望ましい

print(paste0("90%信頼区間は",exp(mean(x0)-mean(x1)+sqrt(1/n0+1/n1)*sigma_hat*qt(0.05,n0+n1-2)),",",exp(mean(x0)-mean(x1)+sqrt(1/n0+1/n1)*sigma_hat*qt(1-0.05,n0+n1-2))))

#推定値は５％信頼区間に含まれるため、棄却されない（有意差がない）

mean(x1)-mean(x0)





```

```{r}

#2薬剤の優劣性の検定(Wilcoxson)　p.70~

data=data.frame(medicine=c("AMPC","S6472"),無効=c(3,8),やや無効=c(8,9),有効=c(30,29),著効=c(22,11),合計=c(63,57))

C=apply(data[,2:(ncol(data)-1)],2,sum)

w=c()

for(j in 1:4){
  
c_sub=C

c_sub[j]=(c_sub[j]+1)/2
  
w=c(w,sum(c_sub[1:j]))  
  
}

n1=data$合計[1];n2=data$合計[2];n=n1+n2

z1=c(data[1,2:5])

W1=sum(w*as.numeric(z1))

E_W1=(n1/n)*sum(w*C)

V_W1=(n1*n2)*(sum(C*w^2)-sum(C*w)^2/n)/(n*(n-1))

sigma_w=sum(C*w^2)/n-(sum(C*w)/n)^2

W=sqrt((n-1)/(n1*n2*sigma_w))*(W1-E_W1)

qnorm(1-0.05/2)

1-pnorm(W)

```

```{r}

#Breslow_Day test

library(dplyr)

data=data.frame(layers=c(1,1,2,2),medicines=c(1,2,1,2),non_effective=c(15,13,7,9),effective=c(21,10,20,34),samples=c(36,23,27,43),effective_rate=c(0.583,0.435,0.741,0.791))

layer1=data %>% filter(layers==1) %>% select(non_effective,effective)

layer2=data %>% filter(layers==2) %>% select(non_effective,effective)

sum=data %>% group_by(layers) %>% summarise(sum=sum(samples))

pthi_MH1=data$non_effective[data$layers==1 & data$medicines==1]*data$effective[data$layers==1 & data$medicines==2]/sum$sum[sum$layers==1]+data$non_effective[data$layers==2 & data$medicines==1]*data$effective[data$layers==2 & data$medicines==2]/sum$sum[sum$layers==2]

pthi_MH2=data$non_effective[data$layers==1 & data$medicines==2]*data$effective[data$layers==1 & data$medicines==1]/sum$sum[sum$layers==1]+data$non_effective[data$layers==2 & data$medicines==2]*data$effective[data$layers==2 & data$medicines==1]/sum$sum[sum$layers==2]

pthi_MH=pthi_MH1/pthi_MH2

#h=1

z12.=sum(layer1[2,]);z1.1=sum(layer1[,1])

z11.=sum(layer1[1,]);

func1=function(x){
  
return((1-pthi_MH)*x^2+((z12.-z1.1)+pthi_MH*(z1.1+z11.))*x-pthi_MH*z11.*z1.1)  
  
}

times=100
newton_raphson=data.frame(times=1:times,value=0,f=0,df=0)
value=50;h=0.01
for(j in 1:nrow(newton_raphson)){
f<-func1(value)
df<-(func1(value+h)-func1(value))/h
  
  newton_raphson$f[j]=f
  newton_raphson$df[j]=df
  value=-f/df+value
  
  newton_raphson$value[j]=value
}

solution=newton_raphson$value[nrow(newton_raphson)]

x1=solution

v1=(1/x1+1/(z11.-x1)+1/(z1.1-x1)+1/(z12.-z1.1+x1))^(-1)

kai_BD1=(layer1[1,1]-x1)^2/v1

#h=2

z12.=sum(layer2[2,]);z1.1=sum(layer2[,1])

z11.=sum(layer2[1,]);

func2=function(x){
  
return((1-pthi_MH)*x^2+((z12.-z1.1)+pthi_MH*(z1.1+z11.))*x-pthi_MH*z11.*z1.1)  
  
}

times=100
newton_raphson=data.frame(times=1:times,value=0,f=0,df=0)
value=50;h=0.01
for(j in 1:nrow(newton_raphson)){
f<-func2(value)
df<-(func2(value+h)-func2(value))/h
  
  newton_raphson$f[j]=f
  newton_raphson$df[j]=df
  value=-f/df+value
  
  newton_raphson$value[j]=value
}

solution=newton_raphson$value[nrow(newton_raphson)]

x2=solution

v2=(1/x2+1/(z11.-x2)+1/(z1.1-x2)+1/(z12.-z1.1+x2))^(-1)


kai_BD2=(layer2[1,1]-x2)^2/v2

kai_BD=kai_BD1+kai_BD2

qchisq(1-0.05,df=1)

#H:p_hij=phi.*ph.j/ph..

x_h_star=c(sum(layer1[,1])*sum(layer1[1,])/as.numeric(sum[1,2]),sum(layer2[,1])*sum(layer2[1,])/as.numeric(sum[2,2]))

z_h11=c(layer1[1,1],layer2[1,1])

n1=sum(apply(layer1,1,sum))

n2=sum(apply(layer2,1,sum))

V_h=c(prod(apply(layer1,1,sum)*apply(layer1,2,sum))/((n1-1)*n1^2),prod(apply(layer2,1,sum)*apply(layer2,2,sum))/((n2-1)*n2^2))

kai_MH=sum(z_h11-x_h_star)^2/sum(V_h)

qchisq(1-0.05,1)

```

```{r}

#3因子交互作用の検定

#H:(a*b*c)=0

m_hij=rep(1,8)

z_hi.=c(apply(layer1,1,sum),apply(layer2,1,sum))

z_h.j=c(apply(layer1,2,sum),apply(layer2,2,sum))

z.ij=c(t(as.matrix(layer1+layer2)))

H=c(1,1,1,1,2,2,2,2)

I=c(1,1,2,2,1,1,2,2)

J=c(1,2,1,2,1,2,1,2)

m_hij_hat=c()

#iterations

times=10

for(k in 1:times){

for(j in 1:length(m_hij)){
  
if(I[j]==1 & H[j]==1){  
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z_hi.[1]/sum(m_hij[1:2])  )
  
}
  
if(I[j]==2 & H[j]==1){  
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z_hi.[2]/sum(m_hij[3:4])  )
  
}  

if(I[j]==1 & H[j]==2){  
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z_hi.[3]/sum(m_hij[5:6]) )
  
}  
  
if(I[j]==2 & H[j]==2){  
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z_hi.[4]/sum(m_hij[7:8]) )
  
}  
  
}  

m_hij=m_hij_hat

m_hij_hat=c()

for(j in 1:length(m_hij)){
  
if(H[j]==1 & J[j]==1){
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z_h.j[1]/sum(m_hij[c(1,3)]) ) 
  
}  
  
if(H[j]==1 & J[j]==2){
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z_h.j[2]/sum(m_hij[c(2,4)]) ) 
  
}  
 
if(H[j]==2 & J[j]==1){
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z_h.j[3]/sum(m_hij[c(5,7)]) ) 
  
}  
  
if(H[j]==2 & J[j]==2){
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z_h.j[4]/sum(m_hij[c(6,8)]) ) 
  
}   
   
}

m_hij=m_hij_hat

m_hij_hat=c()

for(j in 1:length(m_hij)){
  
if(I[j]==1 & J[j]==1){
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z.ij[1]/sum(m_hij[c(1,5)]))  
  
}  
 
if(I[j]==1 & J[j]==2){
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z.ij[2]/sum(m_hij[c(2,6)]))  
  
}  
  
if(I[j]==2 & J[j]==1){
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z.ij[3]/sum(m_hij[c(3,7)]))  
  
}  
 
if(I[j]==2 & J[j]==2){
  
m_hij_hat=c(m_hij_hat,m_hij[j]*z.ij[4]/sum(m_hij[c(4,8)]))  
  
}    
  
   
}

m_hij=m_hij_hat

m_hij_hat=c()

}

z_h=c(c(t(as.matrix(layer1))),c(t(as.matrix(layer2))))

kai=2*sum(z_h*log(z_h)-z_h*log(m_hij))

#yield chisquare distribution for freedom 1

qchisq(1-0.05,1)
  
```