adaptive design methods in clinical trials


#著書：https://www.amazon.co.jp/Adaptive-Methods-Clinical-Chapman-Biostatistics/dp/1584887761/ref=sr_1_1?__mk_ja_JP=%E3%82%AB%E3%82%BF%E3%82%AB%E3%83%8A&keywords=adaptive+design+methods+in+clinical+trials&qid=1581212255&sr=8-1



```{r}

#adaptive design methods in clinical trials

#p.28

#m;the kinds of protocol amendment

#ns;number of patients

ns=30;m=50;m=m-1

x=matrix(sample(0:1,ns*m,replace=T),ncol=ns)

n=prod(dim(x))

mu_j=apply(x,1,mean)

mu_mu=mean(mu_j)

sigma_mu=mean((mu_j-mu_mu)^2)

sigma=sum((x-c(mu_j))^2)/n

#compatible1

sigma_mu*apply(x,1,sum)

sigma*mu_mu

#compatible2

sigma

sigma_mu*ns

x1=x;m1=m;mu_mu1=mu_mu;sigma1=sigma;sigma_mu1=sigma_mu;ns1=ns;n1=n





ns=30;m=5;m=m-1

x=matrix(sample(0:1,ns*m,replace=T),ncol=ns)

n=prod(dim(x))

mu_j=apply(x,1,mean)

mu_mu=mean(mu_j)

sigma_mu=mean((mu_j-mu_mu)^2)

sigma=sum((x-c(mu_j))^2)/n

#compatible1

sigma_mu*apply(x,1,sum)

sigma*mu_mu

#compatible2

sigma

sigma_mu*ns

x2=x;m2=m;mu_mu2=mu_mu;sigma2=sigma;sigma_mu2=sigma_mu;ns2=ns;n2=n


sigma_p=sigma1*sum(m/ns1)/((m1+1)^2)+sigma_mu1/(m1+1)+sigma2*sum(m/ns2)/((m2+1)^2)+sigma_mu2/(m2+1)

#H0:mu_mu1=mu_mu2~N(0,1)

abs(mu_mu1-mu_mu2)/sigma_p

#H0:mu_mu1-mu_mu2<=delta

delta=0.01

(mu_mu1-mu_mu2-delta)/sigma_p

```

```{r}

#2.4 sample size Adjustment

ns=30;m=50;m=m-1

x=matrix(sample(0:1,ns*m,replace=T),ncol=ns)

n=prod(dim(x))

mu_j=apply(x,1,mean)

mu_mu=mean(mu_j)

sigma_mu=mean((mu_j-mu_mu)^2)

sigma=sum((x-c(mu_j))^2)/n


alpha=0.05;beta=0.05;ep=0.1

fun=function(s){
  
s=abs(s)  
  
return((ep-(qnorm(alpha/2)+qnorm(beta))*(sqrt(2/s)*(sigma+s*sigma_mu/(m+1))))^2)
  
}


value=0.1

h=0.01

ite=100

for(j in 1:ite){
  
f<-fun(value)

df<-(fun(value+h)-fun(value))/h

value=value-1/((df)/f)

print(f)  
  
}

value=abs(value)

#R=4*(sigma+sigma_mu)*(qnorm(1-alpha/2)+qnorm(1-beta))^2/(value*ep^2)

R=(1+sigma_mu/sigma)*(1-4*sigma_mu*(qnorm(1-alpha/2)+qnorm(1-beta))^2/((m+1)*ep^2))

```

```{r}

#2.4.2

ns=30;m=50;m=m-1

x=matrix(sample(0:1,ns*m,replace=T),ncol=ns)

n=prod(dim(x))

mu_j=apply(x,1,mean)

mu_mu=mean(mu_j)

sigma_mu=mean((mu_j-mu_mu)^2)

sigma=sum((x-c(mu_j))^2)/n


alpha=0.05;beta=0.05;ep=1;delta=0.5

fun=function(s){
  
s=abs(s)  
  
return((ep-delta-(qnorm(alpha/2)+qnorm(beta))*(sqrt(2/s)*(sigma+s*sigma_mu/(m+1))))^2)
  
}


value=0.1

eta=0.01

h=0.01

ite=100

for(j in 1:ite){
  
f<-fun(value)

df<-(fun(value+h)-fun(value))/h

value=value-1/((df)/f)

print(f)  
  
}

sigma_mu/((ep-delta)^2)

R=1-sigma_mu*(qnorm(1-alpha)+qnorm(1-beta))^2/((m+1)*(ep-delta)^2)

```

```{r}

#2.4.3

ns=30;m=50;m=m-1

x=matrix(sample(0:1,ns*m,replace=T),ncol=ns)

n=prod(dim(x))

mu_j=apply(x,1,mean)

mu_mu=mean(mu_j)

sigma_mu=mean((mu_j-mu_mu)^2)

sigma=sum((x-c(mu_j))^2)/n


alpha=0.05;beta=0.05;ep=1;delta=0.5

fun=function(s){
  
s=abs(s)  
  
return((delta-ep-(qnorm(alpha/2)+qnorm(beta))*(sqrt(2/s)*(sigma+s*sigma_mu/(m+1))))^2)
  
}


value=0.1

eta=0.01

h=0.01

ite=100

for(j in 1:ite){
  
f<-fun(value)

df<-(fun(value+h)-fun(value))/h

value=value-1/((df)/f)

print(f)  
  
}

sigma_mu/((ep-delta)^2)

R=1-sigma_mu*(qnorm(1-alpha)+qnorm(1-beta))^2/((m+1)*(ep-delta)^2)

```

```{r}

#p.38 2.5.1

#m;the kinds of protocol amendment

#ns;number of patients

ns=30;m=50;

x=matrix(sample(0:1,ns*m,replace=T),ncol=ns)

n=prod(dim(x))

mu_j=apply(x,1,mean)

mu_mu=mean(mu_j)

sigma_mu=mean((mu_j-mu_mu)^2)

sigma=sum((x-c(mu_j))^2)/n

#let y_k_bar be the sample mean based on nk clinical data

y_k_bar=c(1)

s_k=c(1)

for(j in 1:(m)){
  
val=rnorm(ns,mu_j[j])

y_k_bar=c(y_k_bar,mean(val))

s_k=c(s_k,var(val))
  
}


X=t(cbind(rep(1,ncol(x)),t(x)))

W=diag(rep(ns,m+1))

beta=solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%y_k_bar

#chisquare-test

s2=sum(ns*s_k)/(sum((m+1)*ns)-(m+1))

qchisq(1-0.05,df=sum((m+1)*ns)-(m+1))

cs=c()

for(j in 1:dim(X)[1]){
  
vec=X[j,]

cs=c(cs,t(vec)%*%solve(t(X)%*%W%*%X)%*%t(t(vec)))
  
}

(c(1,mu_j)-rnorm(length(mu_j)+1,c(1,mu_j)))/as.numeric(sqrt(cs*s2))

```



```{r}

#construction of utility function p.94

library(dplyr)
library(stringr)

data(Indometh)

tau=rnorm(nrow(Indometh))

Indometh=data.frame(Indometh,y=rnorm(nrow(Indometh),0.05))

summary=Indometh%>%group_by(Subject)%>%summarise(n=n())

N=max(summary$n)

m=max(as.numeric(Indometh$Subject))

a_mat=array(0,dim=c(m,N))



for(j in 1:m){

sub=Indometh%>%filter(Subject==j)%>%select(-Subject)

r=ifelse(sub$y>tau[j],1,0)

a=rep(1,nrow(sub));x=sub$conc

fun=function(A){

return(sum(r*log(x*A/sum(x*A))+(1-r)*log(1-x*A/sum(x*A))))

}


h=0.01;eta=0.001

pre_val=-100;val=-97

while(abs(pre_val-val)>0.0001){

vec=a;pre_val=fun(a)

for(i in 1:length(a)){

vec_sub=vec;vec_sub[i]=vec_sub[i]+h

if(vec[i]+eta*(fun(vec_sub)-fun(vec))/h>0){

a[i]=vec[i]+eta*(fun(vec_sub)-fun(vec))/h

}}

val=fun(a)

print(fun(a))

}

a_mat[j,c(1:length(a))]=a

}







```
```{r}

#p.112

library(dplyr)

data(Indometh)

subjects=sort(unique(Indometh$Subject))

lam=c()

for(i in 1:length(subjects)){

sub=Indometh%>%filter(Subject==subjects[i])

times=sub$time;conc=sub$conc

kai2=function(s){

return(sum((conc-s*exp(-s*times))^2/conc))

}


#連続値における適合度検定統計量をもとに推定

ite=1000

eta=0.001

h=0.01

lambda=0.1

val_pre=0.11;val=0.1

while(abs(val_pre-val)>10^(-8))({

val_pre=val

lambda=lambda-eta*(kai2(lambda+h)-kai2(lambda))/h

val=kai2(lambda)

print(kai2(lambda))

})

lam=c(lam,lambda)

}

t0=1;ts=1.1

sigma_vec=lam^2/(1-(exp(lam*t0)-1)/(t0*lam*exp(lam*ts)))


stage=5

#number of patients

Nik=matrix(sample(10:30,length(subjects)*stage,replace=T),ncol=stage)

dik=Nik*c((t0-(exp(lam*t0)-1)/(lam*exp(lam*ts)))/t0)

#2-pairs

delta_mat=array(0,dim=c(length(subjects),length(subjects)))

N_fixed_mat=array(0,dim=c(length(subjects),length(subjects)))

rij=array(0,dim=c(length(subjects),length(subjects)))

alpha=0.05;beta=0.9

for(i in 1:length(subjects)){
for(j in 1:length(subjects)){

delta_mat[i,j]=(lam[i]-lam[j])/sqrt(sigma_vec[i]+sigma_vec[j])

N_fixed_mat[i,j]=(sigma_vec[i]+sigma_vec[j])*(qnorm(1-alpha/2)+qnorm(beta))^2/((lam[i]-lam[j])^2)

}}

#two-samples

N1k=30;N2k=20

d1k=10;d2k=5

r=N2k/N1k

sigma=(N1k+N2k)/(d1k+d2k)

Ik=r*(N1k+N2k)/(sigma*(1+r)^2)

sita=0.1

rnorm(100,sita*Ik,1)

#O'Brien-Fleming

alpha=0.05;s=1;guzai=0.1

a1=2*(1-pnorm(qnorm(alpha)/sqrt(2)))

#Pocock

a2=alpha*log(1+(exp(1)-1)*s)

#Lan-DeMets-Kim

a3=alpha*s^sita

#Hwang-Shih

a4=alpha*(1-exp(guzai*s))/(1-exp(-guzai))


```

```{r}

#p.176

library(dplyr)


n=100

ki=sample(0:1,n,replace=T)

Si=sample(c(c(10:100),Inf),n,replace=T)

Yi=sample(10:100,n,replace=T)

delta=sample(0:1,n,replace=T)



dat=data.frame(ki=ki,Si=Si,Yi=Yi,delta=delta)


log_lik=function(param){

sita=param[1];beta=param[2];et=param[3]

value=0

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

if(dat$Si[j]==Inf){

value=value+dat$delta[j]*(beta*dat$ki[j]+log(dnorm(exp(beta*dat$ki[j]*dat$Yi[j]))))+(1-dat$delta[j])*log(1-pnorm(exp(beta*dat$ki[j]*dat$Yi[j])))

}else{

val1=dnorm(exp(beta*dat$ki[j])*dat$Si[j]+exp(beta*(1-dat$ki[j]))*exp(et*dat$Si[j])*(dat$Yi[j]-dat$Si[j]))

value=value+dat$delta[j]*(beta*(1-dat$ki[j])+log(ifelse(val1>0,(val1),1)))

val2=1-pnorm(exp(beta*dat$ki[j])*dat$Si[j]+exp(beta*(1-dat$ki[j]))*exp(et*dat$Si[j])*(dat$Yi[j]-dat$Si[j]))

value=value+dat$delta[j]*log(ifelse(val2>0,val2,1))

}

}

return(value)

}


ite=1000

eta=0.00001

h=0.01

parameters=rep(0.0001,3)

for(l in 1:ite){

vec=parameters

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

vec_sub=vec;vec_sub[i]=vec_sub[i]+h

parameters[i]=parameters[i]+eta*(log_lik(vec_sub)-log_lik(vec))/h

}

print(log_lik(parameters))

}

```

```{r}

#p.164

n=10;m=5

x=matrix(sample(1:100,n*m,replace=T),ncol=n)

u_hat=apply(x,1,mean)

c=rpois(n,5);c=c-mean(c)

sigma=var(u_hat)

sum(c*u_hat)/sqrt(sigma*sum(c^2/n))



#p.176

n=100

treatment_time=c()

delta=c()

class=5

for(i in 1:class){

time=sample(10:100,1)

treatment_time=c(treatment_time,rpois(n,time))

delta=c(delta,sample(0:1,n,replace=T,prob=c(1/2,1/2)))

}

treat=matrix(treatment_time,ncol=class)

delta=matrix(delta,ncol=class)


log_lik=function(x){

exps=exp(x*treat);exps=t(t(exps)/apply(exps,2,sum))

return(sum(delta*log(exps)))

}


ite=100

eta=0.0001

h=0.01

b=1

for(l in 1:ite){

b=b+eta*(log_lik(b+h)-log_lik(b))/h

print(log_lik(b))

}


#p.177


class=5;n=5

S=sample(10:20,class*n,replace=T)

S_mat=matrix(S,ncol=class)



treatment_time=c()

for(i in 1:class){

time=sample(10:100,1)

treatment_time=c(treatment_time,rpois(n,time))

}

treat=matrix(treatment_time,ncol=class)

delta=matrix(sample(0:1,n*class*class,replace=T),ncol=class)



log_lik=function(x){

beta=x[1];params=x[-1]

params_mat=matrix(params,ncol=class)

w=exp(cbind(S,S^2)%*%params_mat)

exps=w

for(j in 1:class){

exps[(1+n*(j-1)):(n*j),]=exps[(1+n*(j-1)):(n*j),]*exp(beta*treat)

}

exps=log(t(t(exps)/apply(exps,2,sum)))

return(sum(delta*(exps)))

}


eta2=rep(0.01,2*class)

param=c(c(0.01,eta2))

ite=1000000

eta=10^(-11);h=0.01

for(l in 1:ite){

vec=param

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

vec_sub=vec;vec_sub[j]=vec_sub[j]+h

param[j]=param[j]+eta*(log_lik(vec_sub)-log_lik(vec))/h

}

print(log_lik(param))

}

#p.185

library(dplyr)

data(Indometh)

Indometh=Indometh%>%mutate(delta=sample(0:1,nrow(Indometh),replace=T))

n=length(sort(unique(Indometh$Subject)))

t0=min(Indometh$time)

log_lik=function(x){

lam=x

w=c()

for(j in 1:n){

w=c(w,prod(1-lam[j]/lam[-j])^(-1))

}

Fs=c();fs=c()

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

Fs=c(Fs,(t0+sum((exp(-lam*Indometh$time[j])-exp(-lam*(Indometh$time[j]-t0)))*w/lam))/t0)

fs=c(fs,sum(w*(exp(-lam*(Indometh$time[j]-t0))-exp(-lam*Indometh$time[j])))/t0)

}

Ss=1-Fs;

return(sum(Indometh$delta*log(fs))+sum((1-Indometh$delta)*log(Ss)))

}


param=sample(seq(0.01,1,0.01),n)

ite=1000000

eta=10^(-8);h=0.01

for(l in 1:ite){

vec=param

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

vec_sub=vec;vec_sub[j]=vec_sub[j]+h

param[j]=param[j]+eta*(log_lik(vec_sub)-log_lik(vec))/h

}

print(log_lik(param))

}


```