著書：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
# 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
# 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))
# 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)
# 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)
# 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))
# 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
}
# 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))
# 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))
}
# 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))
}