LoginSignup
1
0

More than 3 years have passed since last update.

@kozakai-ryouta(亮太 小酒井)

無作為化比較試験　丹後俊郎先生　オーム社

Last updated at Posted at 2019-06-04
#無作為化比較試験

#同数割り付けの場合の許容区間

n=1000

values=rpois(2,5);p=min(values/sum(values))

samples=sample(c(0:1),n,prob=c(1-p,p),replace = T)

i=500;

average=i/2;var=(i*(n-i))/(4*(n-1))

alpha=0.05

#S_i

upper1=qnorm(1-alpha/2,average,sqrt(var))

under1=qnorm(alpha/2,average,sqrt(var))

#D_ni

upper2=qnorm(1-alpha/2,0,sqrt(4*var))

under2=qnorm(alpha/2,0,sqrt(4*var))


#無作為割り付け

n=1000

values=rpois(2,5);p=min(values/sum(values))

samples=sample(c(0:1),n,prob=c(1-p,p),replace = T)

#D_n:

upper=qnorm(1-alpha/2,0,sqrt(n))

under=qnorm(alpha/2,0,sqrt(n))


#層別無作為化法(Ni=n)

n=100;K=20

N_mat=array(0,dim=c(K,n))

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

values=rpois(2,5);prob=values/sum(values) 

prob=c(1/2,1/2)

N_mat[j,]=sample(c(0:1),n,replace=T,prob=prob)  

}

i=n/2

E_S=i;var_S=(i*(n-i))/(4*(n-1))

E_D=0;var_D=(i*(n-i))/(n-1)

#sum_S

alpha=0.01

value=qnorm(1-0.01,E_S*K,sqrt(var_S*K^2))

apply(N_mat,1,sum)




#母平均差の検定

n1=100;n2=80

X_A=rpois(n1,5);X_B=rpois(n2,9)

s2=((n1-1)*var(X_A)+(n2-1)*var(X_B))/(n1+n2-2)

t=abs(mean(X_A)-mean(X_B))/(sqrt(s2)*sqrt(1/n1+1/n2))

alpha=0.05;beta=0.01;n=100

delta=mean(X_A)-mean(X_B)

fun=function(x){

return(2*((qt(1-alpha/2,2*x-2)+qt(1-beta,2*x-2))/(delta/sqrt(s2)))^2-x  )

}

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

  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)]









#母比率の検定

value1=rpois(2,5);value2=rpois(2,10)

r_A=value1[1];r_B=value2[1]

n_A=sum(value1);n_B=sum(value2)

p_A=r_A/n_A;p_B=r_B/n_B

p_bar=(r_A+r_B)/(n_A+n_B)

t=(p_A-p_B)/((1/n_A+1/n_B)*p_bar*(1-p_bar))

alpha=0.05;beta=0.01

p=(p_A+p_B)/2

n=((qnorm(1-alpha/2)*sqrt(2*p*(1-p))+qnorm(1-beta)*sqrt(p_A*(1-p_A)+p_B*(1-p_B)))/(p_A-p_B))^2






#傾向性検定-量反応関係

K=10

sita=rnorm(K);c_vec=c(1:K)-mean(c(1:K))

r=cor(sita,c_vec)

#平均値の場合

ave_vec=c();var_vec=c();n_vec=c()

for(j in 1:K){

ns=sample(50:100,1)

n_vec=c(n_vec,ns)

y_vec=rnorm(ns,sample(1:10,1),sample(1:10,1))

ave_vec=c(ave_vec,mean(y_vec))  

var_vec=c(var_vec,sum((y_vec-mean(y_vec))^2))

}

N=sum(n_vec)

sigma2=sum(var_vec)/(N-K)

t=sum(c_vec*ave_vec)/sqrt(sigma2*sum(c_vec^2)/(N/K))

alpha=0.05;beta=0.01

n=sum(c_vec^2)*sigma2*(qnorm(1-alpha/2)+qnorm(1-beta))^2/sum(c_vec*ave_vec)^2





#母比率の場合

K=10

r_vec=c();n_vec=c()

for(j in 1:K){

value=rpois(2,sample(1:10,1))  

r_vec=c(r_vec,value[1])

n_vec=c(n_vec,sum(value))


}

p=r_vec/n_vec

n=sum(n_vec)/K

p_bar=sum(r_vec)/sum(n_vec)

t=sum(c_vec*p)/sqrt(sum(c_vec^2)*p_bar*(1-p_bar)/n)

alpha=0.05;beta=0.01

n=((qnorm(1-alpha/2)*sqrt((sum(p)/K)*(1-sum(p)/K)*sum(c_vec^2))+qnorm(1-beta)*sqrt(sum(p*(1-p)*c_vec^2)))/sum(c_vec*p))^2





#log-rank test


#compared-hazard model

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

param=cov(data$num,log(data$y))/var(data$num,data$num)

intercept=mean(log(data$y))-param*mean(data$num)

predict=param*data$num+intercept

y=log(data$y)-predict

X=as.matrix(data[,colnames(data) %in% c("x1","x2")])

param2=solve(t(X)%*%X)%*%t(X)%*%y

predict_data=data.frame(num=1:nrow(data),y=log(data$y),predict=predict+X%*%param2)

plot(predict_data$num,(predict_data$y),col=2,type="o",ylim=c(min(log(data$y),predict_data$predict),max(log(data$y),predict_data$predict)),xlab="index",ylab="values")

par(new=T)

plot(predict_data$num,predict_data$predict,col=3,type="o",ylim=c(min(log(data$y),predict_data$predict),max(log(data$y),predict_data$predict)),xlab="index",ylab="values")



#log-rank

n=100

data1=data.frame(t=1:n,x=0,event=sample(1:10,n,replace = T),survive=sample(1:10,n,replace = T)) %>% mutate(R=event+survive)

data2=data.frame(t=1:n,x=1,event=sample(1:10,n,replace = T),survive=sample(1:10,n,replace = T)) %>% mutate(R=event+survive)

data=rbind(data1,data2)

n_data=data %>% group_by(t) %>% summarise(n=sum(R))

E_H0_d1j=c();var_H0_d1j=c();d1j=c();r_vec=c()

for(j in 1:n){

n1j=data2$R[data2$t==j];n0j=data1$R[data1$t==j];dj=data2$event[data2$t==j]+data1$event[data1$t==j];nj=n1j+n0j  

E_H0_d1j=c(E_H0_d1j,n1j*dj/(n1j+n0j))  

var_H0_d1j=c(var_H0_d1j,((nj-n1j)/(nj-1))*n1j*(dj/nj)*(1-dj/nj)) 

d1j=c(d1j,data2$event[data2$t==j])

r_vec=c(r_vec,n1j/n0j)

}

kai2=sum(d1j-E_H0_d1j)^2/sum(var_H0_d1j)

qchisq(1-0.05,1)

#統計量

t=sum(d1j-E_H0_d1j)

E_H0_t=0

var_H0_t=sum(r_vec/(r_vec+1)^2)

#平均生存率(x=0)

plot(data1$survive/data1$R)

#平均生存率(x=1)

plot(data2$survive/data2$R)

sita_hat=mean(data2$survive/data2$R)/mean(data1$survive/data1$R)

alpha=0.05;beta=0.01

e=(qnorm(1-alpha/2)+qnorm(1-beta))*((sita_hat+1)/(sita_hat-1))^2

N=e/(2-mean(data2$survive/data2$R)-mean(data1$survive/data1$R))





#非劣性の検定(Dunneett-Gent(1977))

delta=0.1;alpha=0.05;beta=0.01

value1=rpois(2,5);value2=rpois(2,9)

r_A=value1[1];r_B=value2[1]

n_A=sum(value1);n_B=sum(value2)

p_A=r_A/n_A;p_B=r_B/n_B

p_B2=(r_A+r_B+n_A*delta)/(n_A+n_B)

S.E.=sqrt((p_B2-delta)*(1-p_B2+delta)/n_A+(p_B2*(1-p_B2))/n_B)

#正規分布に従うN(0,1)

Z=(p_A-p_B)/S.E.


#H0:p_A=p_B+delta

a=n_A+n_B

b=-(n_B+n_A+r_B+r_A+delta*(n_A+2*n_B))

c=n_B*delta^2+delta*(2*r_B+n_A+n_B)+r_B+r_A

d=-r_B*delta*(1+delta)

v=b^3/(27*a^3)-b*c/(6*a^2)+d/(2*a)

u=sign(v)*sqrt(b^2/(9*a^2)-c/(3*a))

w=(pi+acos(v/(u^3)))/3

p_B_star=2*u*cos(w)-b/(3*a)

p_B_bar=p_B+(p_A-p_B+delta)/2

#sample size

N=((qnorm(1-alpha)*sqrt((p_B_bar-delta)*(1-p_B_bar+delta)+p_B_bar*(1-p_B_bar))+qnorm(1-beta)*sqrt(p_A*(1-p_A)+p_B*(1-p_B)))/(p_A-p_B+delta))^2

#p_B_starを利用した場合のサンプルサイズ

a=2;

b=-2*p_B-2-3*delta-(p_A-p_B)

c=delta^2+2*(1+p_B)*delta+2*p_B+(p_A-p_B)

d=-p_B*delta*(1+delta)

v=b^3/(27*a^3)-b*c/(6*a^2)+d/(2*a)

u=sign(v)*sqrt(b^2/(9*a^2)-c/(3*a))

w=(pi+acos(v/(u^3)))/3

p_B_star2=2*u*cos(w)-b/(3*a)

#sample size

N2=((qnorm(1-alpha)*sqrt((p_B_star2-delta)*(1-p_B_star2+delta)+p_B_star2*(1-p_B_star2))+qnorm(1-beta)*sqrt(p_A*(1-p_A)+p_B*(1-p_B)))/(p_A-p_B+delta))^2



#同等性の検定

n=10

delta=sample(seq(0.01,0.1,0.01),n)

p_A_vec=c();p_B_vec=c();r_A_vec=c();r_B_vec=c();n_A_vec=c();n_B_vec=c()

for(j in 1:n){

value1=rpois(2,5);value2=rpois(2,9)

r_A=value1[1];r_B=value2[1]

n_A=sum(value1);n_B=sum(value2)

p_A=r_A/n_A;p_B=r_B/n_B 

p_A_vec=c(p_A_vec,p_A);p_B_vec=c(p_B_vec,p_B)

r_A_vec=c(r_A_vec,r_A);r_B_vec=c(r_B_vec,r_B)

n_A_vec=c(n_A_vec,n_A);n_B_vec=c(n_B_vec,n_B)

}

#H0:pthi=1;H1:pthi>1

pthi=p_A_vec*(1-p_B_vec+delta)/((1-p_A_vec)*(p_B_vec-delta))

p_B_star_vec=c()

for(j in 1:n){

fun=function(x){

z<-log((((pthi[j]*(x-delta[j]))/(pthi[j]*(x-delta[j])+1-x+delta[j]))^r_A_vec[j])*(((1-x+delta[j])/(pthi[j]*(x-delta[j])+1-x+delta[j]))^(n_A_vec[j]-r_A_vec[j]))*((x^r_B_vec[j])*(1-x)^(n_B_vec[j]-r_B_vec[j])) ) 

return(z)  

}

value=0.001

newton_raphson=data.frame(x=seq(0,1,value),value=fun(seq(0,1,value))) %>% filter(is.nan(value)!=T)

solution=newton_raphson$x[newton_raphson$value==max(newton_raphson$value)]

p_B_star_vec=c(p_B_star_vec,min(solution))

}

U=sum(r_A_vec-n_A_vec*(p_B_star_vec-delta))

var_H0=sum((n_A_vec*n_B_vec*(p_B_star_vec-delta)^2)*((1-p_B_star_vec+delta)^2)/(n_A_vec*p_B_star_vec*(1-p_B_star_vec)+n_B_vec*(p_B_star_vec-delta)*(1-p_B_star_vec+delta)))

Z=U/sqrt(var_H0)






library(dplyr)

#非劣性の検定(H0:p_A<=p_B-delta)

mat=matrix(c(3,5,7,8),ncol=2)

colnames(mat)=c("対照薬（有効）","対照薬（無効）")

rownames(mat)=c("対照薬（有効）","対照薬（無効）")

a=mat[1,1];b=mat[1,2];c=mat[2,1];d=mat[2,2];n=sum(mat)

q12=b/n;q21=c/n;delta=0.1

t=(b-c+n*delta)/n

var_t=(q12+q21-delta^2)/n


B=(-b-c-(2*n-b+c)*delta)/n

C=(c*delta*(delta+1))/n

q21_hat=(sqrt(B^2-8*C)-B)/4

#漸近正規N(0,1)

Z=(b-c+n*delta)/sqrt(n*(2*q21_hat-delta-delta^2))


#sample size

lambda=apply(mat,1,sum)[1]/n-apply(mat,2,sum)[1]/n

B_star=-lambda*(1-delta)-2*(q21+delta)

C_star=q21*delta*(delta+1)

q21_star=(sqrt(B_star^2-8*C_star)-B_star)/4

alpha=0.05;beta=0.01

N=((2*q21_star-delta-delta^2)*(qnorm(1-alpha)/sqrt(2*q21_star-delta-delta^2)+qnorm(1-beta)/sqrt(2*q21+lambda-lambda^2)))^2



#Pocok method 1

library(dplyr)

times=10;n=100

dat_A=data.frame(t=1:times,X=0)

dat_B=data.frame(t=1:times,X=0)

S=c()

for(j in 1:times){

X_A=(rnorm(n,sample(1:10,1),sample(1:10,1))) 

X_B=(rnorm(n,sample(1:10,1),sample(1:10,1)))  

dat_A$X[j]=mean(X_A) 

dat_B$X[j]=mean(X_B) 


S=c(S,((n-1)*var(X_A)+(n-1)*var(X_B))/(n+n-2))

}

dat=data.frame(times=1:times,X_A=dat_A$X,X_B=dat_B$X,S=S)

t=c()

for(j in 1:times){

t=c(t,(dat_A$X[j]-dat_B$X[j])/(S[j]*sqrt(2/n))) 

}

a=sort(sample(1:10,times,replace = T))

alpha=1-(pnorm(a[times],0,times)-pnorm(-a[times],0,times))

prob=c()

pthi=2

for(j in 2:times){

prob=c(prob,1-(1-(pnorm(a[j],pthi*j,j)-pnorm(-a[j],pthi*j,j)) ) )

}

ASN=2*n*(1+sum(prob))




#Pocok method 2

library(dplyr)

times=10;n=100

dat_A=data.frame(t=1:times,X=0)

dat_B=data.frame(t=1:times,X=0)

S=c();pthi=c()

for(j in 1:times){

val_A=rpois(2,sample(1:10,1));val_B=rpois(2,sample(1:10,1))

X_A=val_A[1]/sum(val_A) 

X_B=val_B[1]/sum(val_B)

dat_A$X[j]=X_A

dat_B$X[j]=X_B 

X=(X_A+X_B)/2

S=c(S,sqrt(2*X*(1-X)/n))

}

dat=data.frame(times=1:times,X_A=dat_A$X,X_B=dat_B$X,S=S)

t=c()

for(j in 1:times){

t=c(t,(dat_A$X[j]-dat_B$X[j])/(S[j])) 

}

a=sort(sample(1:10,times,replace = T))

alpha=1-(pnorm(a[times],0,times)-pnorm(-a[times],0,times))

prob=c()

pthi=2

for(j in 2:times){

prob=c(prob,1-(1-(pnorm(a[j],pthi*j,j)-pnorm(-a[j],pthi*j,j)) ) )

}

ASN=2*n*(1+sum(prob))




#alpha消費関数

alpha=0.05

a_t=function(t){

return(ifelse(t>0,2-2*pnorm(qnorm(1-alpha/2)/sqrt(t)),0))  

}

a_t_data=data.frame(t=seq(0,1,0.01),a_t=a_t(seq(0,1,0.01))) %>% mutate(c_t=sqrt(t)*qnorm(1-a_t))

plot(a_t_data$t,a_t_data$a_t,type="p",col=2)





#alpha消費関数に基づく逐次検定統計量の同時分布

library(dplyr)

times=10;n_vec_A=c();n_vec_B=c()

dat_A=data.frame(t=(1:times)/times,X=0)

dat_B=data.frame(t=(1:times)/times,X=0)

S=c()

X_A_vec=c();X_B_vec=c()

sigma_A=c();sigma_B=c()

for(j in 1:times){

N_A=sample(10:100,1)

N_B=sample(10:100,1)

n_vec_A=c(n_vec_A,N_A);n_vec_B=c(n_vec_B,N_B)

X_A=(rnorm(N_A,sample(1,1),sample(1,1))) 

X_B=(rnorm(N_B,sample(1,1),sample(1,1)))  

X_A_vec=c(X_A_vec,X_A)

X_B_vec=c(X_B_vec,X_B)

sigma_A=c(sigma_A,var(X_A_vec))

sigma_B=c(sigma_B,var(X_B_vec))

dat_A$X[j]=sum(X_A) 

dat_B$X[j]=sum(X_B) 

S=c(S,((N_A-1)*var(X_A)+(N_B-1)*var(X_B))/(N_A+N_B-2))

}

dat_A$X=dat_A$X/cumsum(n_vec_A)

dat_B$X=dat_B$X/cumsum(n_vec_B)

dat=data.frame(times=(1:times)/times,W_A=dat_A$X,W_B=dat_B$X,S=(dat_A$X-dat_B$X)/sqrt(S),v=(1/cumsum(n_vec_A)+1/cumsum(n_vec_B))) %>% mutate(U=S/sqrt(1/cumsum(n_vec_A)+1/cumsum(n_vec_B)))


#alpha消費関数

alpha=0.05

a_t=function(t){

return(ifelse(t>0,2-2*pnorm(qnorm(1-alpha/2)/sqrt(t)),0))  

}

a_t_data=data.frame(t=seq(0,1,0.01),a_t=a_t(seq(0,1,0.01))) %>% mutate(c_t=sqrt(t)*qnorm(1-a_t))

plot(a_t_data$t,a_t_data$a_t,type="p",col=2)



library(mvtnorm)

t1=0.3;t2=0.9

b_vec=c();alpha_vec=c()

b1=qnorm(1-a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==t1]/2)

b_vec=c(b_vec,b1);alpha_vec=c(alpha_vec,a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==t1]/2)

t_s=seq(t1,t2,0.1)

for(j in 2:(length(t_s))){

t_vec=t(t(t_s[1:j]))

cov_U=sqrt(t_vec%*%t(1/t_vec))

cov_U[lower.tri(cov_U)]=0

cov_U=cov_U+t(cov_U);diag(cov_U)=diag(cov_U)/2

b_vec=c(b_vec,qmvnorm( a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==floor(t_s[j]*10^2)/(10^2)]- a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==floor(t_s[j-1]*10^2)/(10^2)],mean = rep(0,ncol(cov_U)),sigma=cov_U)$quantile)

alpha_vec=c(alpha_vec,a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==floor(t_s[j]*10^2)/(10^2)]-a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==floor(t_s[j-1]*10^2)/(10^2)])

}

b_vec=abs(b_vec)

W_A=dat$W_A[2:times];W_B=dat$W_B[2:times]

sigma_hat=sqrt((sigma_A*(cumsum(n_vec_A)-1)+sigma_B*(cumsum(n_vec_B)-1))/(cumsum(n_vec_A)+cumsum(n_vec_B)-2))

under=W_A-W_B-b_vec*sqrt(sigma_hat[2:times])*sqrt(dat$v[2:times])

upper=W_A-W_B+b_vec*sqrt(sigma_hat[2:times])*sqrt(dat$v[2:times])

#6.7.1　反復信頼区間

plot((2:times)/times,under,type="o",col=2,ylim=c(min(under),max(upper)),xlab="times",ylab="values")

par(new=T)

plot((2:times)/times,upper,type="o",col=3,ylim=c(min(under),max(upper)),xlab="times",ylab="values")

res=data.frame(t=t_s,b_vec,alpha_vec,p_val=0)

M=3000

for(l in 1:M){

n_vec_A=c();n_vec_B=c()

dat_A=data.frame(t=(1:times)/times,X=0)

dat_B=data.frame(t=(1:times)/times,X=0)

S=c()

X_A_vec=c();X_B_vec=c()

sigma_A=c();sigma_B=c()

for(j in 1:times){

N_A=sample(10:100,1)

N_B=sample(10:100,1)

n_vec_A=c(n_vec_A,N_A);n_vec_B=c(n_vec_B,N_B)

X_A=(rnorm(N_A,sample(1,1),sample(1,1))) 

X_B=(rnorm(N_B,sample(1,1),sample(1,1)))  

X_A_vec=c(X_A_vec,X_A)

X_B_vec=c(X_B_vec,X_B)

sigma_A=c(sigma_A,var(X_A_vec))

sigma_B=c(sigma_B,var(X_B_vec))

dat_A$X[j]=sum(X_A) 

dat_B$X[j]=sum(X_B) 

S=c(S,((N_A-1)*var(X_A)+(N_B-1)*var(X_B))/(N_A+N_B-2))

}

dat_A$X=dat_A$X/cumsum(n_vec_A)

dat_B$X=dat_B$X/cumsum(n_vec_B)

dat=data.frame(times=(1:times)/times,W_A=dat_A$X,W_B=dat_B$X,S=(dat_A$X-dat_B$X)/sqrt(S),v=(1/cumsum(n_vec_A)+1/cumsum(n_vec_B))) %>% mutate(U=S/sqrt(1/cumsum(n_vec_A)+1/cumsum(n_vec_B)))

vecs=c()

for(k in 1:length(t_s)){

 u_vec=abs(dat$U[as.character(floor(dat$times*10^2)/(10^2)) %in% as.character(floor(t_s*10^2)/(10^2))])[1:k]

 vecs=c(vecs,ifelse(sum(u_vec[1:(k-1)]<=b_vec[1:(k-1)])==(k-1) & u_vec[k]>b_vec[k],1,0))

}

res$p_val=res$p_val+vecs

print(res$p_val/rep(l,nrow(res)))

}

res$p_val=res$p_val/M

#2-stage design

#fisher p-value

-2*sum(log(res$p_val))>qchisq(1-0.01,2*nrow(res))

#inverse normal

sum(qnorm(1-res$p_val))/sqrt(nrow(res))>qnorm(1-0.01)



#sample size

times=10

n_vec_A=c();n_vec_B=c()

dat_A=data.frame(t=(1:times)/times,X=0)

dat_B=data.frame(t=(1:times)/times,X=0)

S=c()

X_A_vec=c();X_B_vec=c()

sigma_A=c();sigma_B=c()

for(j in 1:times){

N_A=sample(10:100,1)

N_B=sample(10:100,1)

n_vec_A=c(n_vec_A,N_A);n_vec_B=c(n_vec_B,N_B)

X_A=(rnorm(N_A,sample(1,1),sample(1,1))) 

X_B=(rnorm(N_B,sample(1,1),sample(1,1)))  

X_A_vec=c(X_A_vec,X_A)

X_B_vec=c(X_B_vec,X_B)

sigma_A=c(sigma_A,var(X_A_vec))

sigma_B=c(sigma_B,var(X_B_vec))

dat_A$X[j]=sum(X_A) 

dat_B$X[j]=sum(X_B) 

S=c(S,((N_A-1)*var(X_A)+(N_B-1)*var(X_B))/(N_A+N_B-2))

}

dat_A$X=dat_A$X/cumsum(n_vec_A)

dat_B$X=dat_B$X/cumsum(n_vec_B)

dat=data.frame(times=(1:times)/times,W_A=dat_A$X,W_B=dat_B$X,S=(dat_A$X-dat_B$X)/sqrt(S),v=(1/cumsum(n_vec_A)+1/cumsum(n_vec_B))) %>% mutate(U=S/sqrt(1/cumsum(n_vec_A)+1/cumsum(n_vec_B)))



w=sqrt(c(dat$times[1],dat$times[2:nrow(dat)]-dat$times[1:(nrow(dat)-1)]))



sigma_hat=sqrt(((sum(n_vec_A)-1)*var(X_A_vec)+(sum(n_vec_B)-1)*var(X_B_vec))/(sum(n_vec_A)+sum(n_vec_B)-2))

t_vec=(dat_A$X-dat_B$X)/(sigma_hat*sqrt(1/n_vec_A+1/n_vec_B))

#alpha消費関数

alpha=0.05

a_t=function(t){

return(ifelse(t>0,2-2*pnorm(qnorm(1-alpha/2)/sqrt(t)),0))  

}

a_t_data=data.frame(t=seq(0,1,0.01),a_t=a_t(seq(0,1,0.01))) %>% mutate(c_t=sqrt(t)*qnorm(1-a_t))

plot(a_t_data$t,a_t_data$a_t,type="p",col=2)



library(mvtnorm)

t1=0.1;t2=1

b_vec=c();alpha_vec=c()

b1=qnorm(1-a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==t1]/2)

b_vec=c(b_vec,b1);alpha_vec=c(alpha_vec,a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==t1]/2)

t_s=seq(t1,t2,0.1)

for(j in 2:(length(t_s))){

t_vec=t(t(t_s[1:j]))

cov_U=sqrt(t_vec%*%t(1/t_vec))

cov_U[lower.tri(cov_U)]=0

cov_U=cov_U+t(cov_U);diag(cov_U)=diag(cov_U)/2

b_vec=c(b_vec,qmvnorm( a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==floor(t_s[j]*10^2)/(10^2)]- a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==floor(t_s[j-1]*10^2)/(10^2)],mean = rep(0,ncol(cov_U)),sigma=cov_U)$quantile)

alpha_vec=c(alpha_vec,a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==floor(t_s[j]*10^2)/(10^2)]-a_t_data$a_t[floor(a_t_data$t*10^2)/(10^2)==floor(t_s[j-1]*10^2)/(10^2)])

}

b_vec=abs(b_vec)

a_vec=sqrt(c(1:nrow(dat)))*b_vec

delta=0.1;beta=0.01

n_vec=2*((w*qnorm(1-beta)+a_vec*sqrt(dat$times/c(1:nrow(dat)))-c(0,cumsum(w*t_vec)[1:(length(w)-1)]))/(w*(delta/sigma_hat)))^2



#cross over test
1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0