医学への統計学　丹後俊郎先生　朝倉出版

```{r}

#仮説検定の考え方　p.44

data=data.frame(No=c(1:10),原因不明の急性出血性腸炎患者=c(8.4,9.4,7.7,9.4,7.8,8.6,8.9,8.2,9.6,8.8),薬剤性腸炎患者=c(7.3,8.1,7.7,6.4,6.8,7.2,7.8,9,NA,NA))

A=data$原因不明の急性出血性腸炎患者
B=data$薬剤性腸炎患者[is.na(data$薬剤性腸炎患者)!=T]

statistics=array(0,dim=c(4,2))

statistics[1,1]=mean(A);statistics[1,2]=mean(B)
statistics[2,1]=sd(A);statistics[2,2]=sd(B)
statistics[3,1]=var(A);statistics[3,2]=var(B)
statistics[4,]=c(10,8)

rownames(statistics)=c("mean","sd","var","numbers")

colnames(statistics)=c("原因不明の急性出血性腸炎患者","薬剤性腸炎患者")

delta=statistics[1,1]-statistics[1,2]

SE_delta=sqrt((1/statistics[4,1]+1/statistics[4,2])*((statistics[4,1]-1)*statistics[3,1]+(statistics[4,2]-1)*statistics[3,2])/(statistics[4,1]+statistics[4,2]-2))

t=delta/SE_delta

alpha=0.05

qt(alpha/2,df=sum(statistics[4,])-1)


```


```{r}



data=data.frame(dead=c(0:7),observations=c(31,44,33,8,3,1,0,1),expect=c(32.8,42.8,28,12.2,4,1,0.2,0))

A_day_of_expect_died_people=sum(data$dead*data$observations)/sum(data$observations)

lambda=A_day_of_expect_died_people

data=data %>% mutate(predict_ave=sum(data$observations)*(lambda^2)*exp(-lambda)/gamma(dead+1))

#negative binomial

lambda_hat=mean(data$observations)

pthi_hat=lambda^2/(var(data$observations)-lambda)

pthi_exam=pthi_hat # The difference of individuals

negative_binom=function(r,pthi,ave_observations){
  
z<-((pthi/(pthi+ave_observations))^pthi)*((ave_observations/(pthi+ave_observations))^r)*prod((c(0:(r-1))+pthi)/c(1:r)  )

return(z)

}

#average

E=0

for(j in 1:1000){
  
E=E+j*negative_binom(j,pthi_exam,lambda)  
  
}

#Check result of expectation value

floor(E*1000)/1000==floor(lambda*1000)/1000

#variance

V=0

for(j in 1:10000){
  
V=V+((j-E)^2)*negative_binom(j,pthi_exam,lambda)  
  
}



#Check result of standard deviation

sqrt(V)

sqrt(E+(E^2)/pthi_exam)

#MLE

samples=sample(c(1:10),50,replace=T)

sample_summary=data.frame(sample=samples) %>% group_by(sample) %>% summarise(n=n())

lambda=mean(samples)

pois=function(x,lambda){
  
z<-(lambda^x)*exp(-lambda)/(gamma(x+1))

return(z)
  
}

value=0

for(j in 1:100){

value=value+((pois(j+h,lambda)-pois(j,lambda))/h)^2

}

n=length(samples)

var_estimate=1/(n*value)

#parameter:lambda (confidence interval)

upper=lambda+2*sqrt(var_estimate)

under=lambda-2*sqrt(var_estimate)

dnorm(samples,mean=lambda,sd=sqrt(var_estimate))

#t-test;f-test

sample_A=sample(c(1:10),10,replace=T)

sample_B=sample(c(1:10),5,replace=T)

SA=var(sample_A);SB=var(sample_B)

f=SA/SB;

qf(0.025,length(sample_A)-1,length(sample_B)-1)

if(f<qf(0.025,length(sample_A)-1,length(sample_B)-1)){

t=(mean(sample_A)-mean(sample_B))/sqrt((1/length(sample_A)+1/length(sample_B))*(((length(sample_A)-1)*var(sample_A)+(length(sample_B)-1)*var(sample_B))/(length(sample_A)+length(sample_B)-2)))

print(paste0("test;",t,": statics;",qt(0.025,df=length(sample_A)+length(sample_B)-2)))

}else{
  
 t=(mean(sample_A)-mean(sample_B))/sqrt(SA/length(sample_A)+SB/length(sample_B))
 
 c=(SA^2/length(sample_A))/((SA^2/length(sample_A))+(SB^2/length(sample_B)))
 
 v=1/((c^2)/(length(sample_A)-1)+((1-c)^2)/(length(sample_B)-1))
 
 print(paste0("test;",t,": statistics;",qt(0.025,df=floor(v))))
  
}


#Grubbs-Smirnov test

data=data.frame(no=c(1:10),T_Chol=c(120,125,118,114,130,154,124,122,117,126))

X_bar=mean(data$T_Chol)

S=sd(data$T_Chol)

t=(data$T_Chol-X_bar)/S



```

```{r}
#wilcoxson test1



library(dplyr)

n=100

data=data.frame(num=1:n,values=rnorm(n,mean=0,sd=2))

t=0

for(j in 1:nrow(data)){
  
value=data$values[j]

t=t+length(data$values[abs(data$values)<abs(value)])*sign(value)
  
  
}

z=sqrt((6/(n*(n+1)*(2*n+1))))

abs(qnorm(0.05/2))

#plot

plot(c(1:length(data$values[data$values<0])),(abs(data$values[data$values<0])),xlab="num",ylab="values",type="p",col=2)



plot(c(1:length(data$values[data$values>0])),(abs(data$values[data$values>0])),xlab="num",ylab="values",type="p",col=3)

data1=data.frame(num=c(1:length(data$values[data$values<0])),values=(abs(data$values[data$values<0])))

data2=data.frame(num=c(1:length(data$values[data$values>0])),values=(abs(data$values[data$values>0])))

data1=data1[order(data1$values),]

data2=data2[order(data2$values),]

```

```{r}
#wilcoxson test2



library(dplyr)

n=50

data1=data.frame(num=1:n,values=sample(1:5,n,replace=T,prob=c(0.0,0.1,0.1,0.1,0.7)))

data2=data.frame(num=1:n,values=-sample(1:5,n,replace=T,prob=c(0.3,0.2,0.1,0.2,0.2)))

data=rbind(data1,data2)

t=0

for(j in 1:nrow(data)){
  
value=data$values[j]

t=t+length(data$values[abs(data$values)<abs(value)])*sign(value)
  
  
}



n=nrow(data)

z=sqrt((6/(n*(n+1)*(2*n+1))))*t

abs(qnorm(0.05/2))

#plot

plot(c(1:length(data$values[data$values<0])),(abs(data$values[data$values<0])),xlab="num",ylab="values",type="p",col=2)



plot(c(1:length(data$values[data$values>0])),(abs(data$values[data$values>0])),xlab="num",ylab="values",type="p",col=3)

data1=data.frame(num=c(1:length(data$values[data$values<0])),values=(abs(data$values[data$values<0])))

data2=data.frame(num=c(1:length(data$values[data$values>0])),values=(abs(data$values[data$values>0])))

data1=data1[order(data1$values),]

data2=data2[order(data2$values),]

z

```

```{r}
#Grubbs Smirnov test

#H0:飛び離れたデータは存在しない

data=data.frame(num=c(1:10),values=c(120,125,118,114,130,154,124,122,117,126))

average=mean(data$values)

t=(max(data$values)-average)/sd(data$values)

#Grubbs smirnov test value(0.01):2.41

if(t>2.41){
 
print("Hypothesis is rejected")   
  
}

```

```{r}

#wilcoxsonによる分布の対称性と分布の類似性の検定

patient=c(10,13,7,16,12,12,10,10,4,17,14,6,11,12,7,9,12,9)

able_bodies_person=c(7,8,8,6,7,11,10,8,9)

patient_data=data.frame(num=c(1:length(patient)),values=patient)

patient_data=patient_data %>% group_by(values) %>% summarise(n=n())

able_bodies_person_data=data.frame(num=c(1:length(able_bodies_person)),values=as.numeric(able_bodies_person))

able_bodies_person_data=able_bodies_person_data %>% group_by(values) %>% summarise(n=n())


plot(patient_data$values,patient_data$n,type="l")

plot(able_bodies_person_data$values,able_bodies_person_data$n,type="l")

patient_ave=mean(patient)

patient_var=var(patient)

able_bodies_person_ave=mean(able_bodies_person)

able_bodies_person_var=var(able_bodies_person)

#等分散性の検定（棄却）

f=patient_var/able_bodies_person_var

f>qf(0.025,length(patient)-1,length(able_bodies_person)-1)

t=abs(able_bodies_person_ave-patient_ave)/sqrt(patient_var/length(patient)+able_bodies_person_var/length(able_bodies_person))

#平均値の検定（棄却）

t>qt(0.025,df=length(patient)+length(able_bodies_person))

#wilcoxson test(対称性)

data=data.frame(num=c(1:(length(patient)+length(able_bodies_person))),values=sort(c(patient,able_bodies_person)))

num_order_data=data %>% group_by(values) %>% summarise(order=mean(num))

data=left_join(data,num_order_data)

able_bodies_person_order=c()

for(j in 1:length(able_bodies_person)){
  
able_bodies_person_order=c(able_bodies_person_order,num_order_data$order[num_order_data$values==able_bodies_person[j]])
  
}

t=sum(able_bodies_person_order)

#正規分布での検定

m=length(patient)+length(able_bodies_person)

mean_t=(1/4)*m*(m+1)

var_t=(1/24)*m*(m+1)*(2*m+1)

(t-mean_t+1/2)/sqrt(var_t)

#この結果より、patientの平均の方が大きい

```

```{r}
#p.94

library(dplyr)

#千葉県のデータ：X;NO2濃度、Y;有症率

data1=data.frame(area_no=c(1:11),X=c(0.025,0.03,0.016,0.028,0.021,0.041,0.012,0.015,0.03,0.018,0.023),Y=c(4.1,6.8,2.6,11.5,2.8,10.7,2.7,2.8,3.5,5,3.6))

SSXY1=sum(data1$X*data1$Y)-sum(data1$X)*sum(data1$Y)/nrow(data1)

SSX1=sum(data1$X^2)-sum(data1$X)^2/nrow(data1)

SSY1=sum(data1$Y^2)-sum(data1$Y)^2/nrow(data1)

X1_mean=mean(data1$X)

Y1_mean=mean(data1$Y)

b1=SSXY1/SSX1

a1=Y1_mean-X1_mean*b1

r1=cor(data1$X,data1$Y)

#岡山県のデータ：X;NO2濃度、Y;有症率

data2=data.frame(area_no=c(1:10),X=c(0.024,0.033,0.03,0.017,0.022,0.027,0.016,0.017,0.02,0.022),Y=c(6.8,8.7,11.4,4.5,5.6,10.5,5.2,6.8,6,8))

SSXY2=sum(data2$X*data2$Y)-sum(data2$X)*sum(data2$Y)/nrow(data2)

SSX2=sum(data2$X^2)-sum(data2$X)^2/nrow(data2)

SSY2=sum(data2$Y^2)-sum(data2$Y)^2/nrow(data2)

X2_mean=mean(data2$X)

Y2_mean=mean(data2$Y)

b2=SSXY2/SSX2

a2=Y2_mean-X2_mean*b2

r2=cor(data2$X,data2$Y)

#H0:b1=b2

N=nrow(data1)+nrow(data2)

delta1=SSY1-(SSXY1^2)/SSX1+SSY2-(SSXY2^2)/SSX2

delta2=SSY1+SSY2-((SSXY1+SSXY2)^2)/(SSX1+SSX2)

f1=(delta2/delta1-1)*((N-4))

if(f1<qf(1-0.05,1,N-4)){
  
  print("棄却されない")
  
}else{
  
  print("棄却")
  
}

#K0:a1=a2

data=rbind(data1,data2)

SSXY=sum(data$X*data$Y)-sum(data$X)*sum(data$Y)/nrow(data)

SSX=sum(data$X^2)-sum(data$X)^2/nrow(data)

SSY=sum(data$Y^2)-sum(data$Y)^2/nrow(data)

delta3=SSY-SSXY^2/SSX

f2=(delta3/delta2-1)*(nrow(data)-3)

if(f2<qf(1-0.01,1,N-3)){
  
  print("棄却されない")
  
}else{
  
  print("棄却")
  
}

b_E=(SSXY1+SSXY2)/(SSX1+SSX2)

a1=Y1_mean-b_E*X1_mean

a2=Y2_mean-b_E*X2_mean

```