###############################################################################
#rm(list=ls())
#data=read.table("covid19.csv", sep=",", head=T)
data = read_excel("C:/Users/mindy/Dropbox/Novel Corona Virus/Multi-center Data Analysis/multi_center.xlsx")
data = data %>% filter(!(is.na(Time2Death)&is.na(Time2Discharge)))
death=data$Death
arbidol=data$Arbidol
kaletra=data$Kaletra
tamiflu=data$Oseltamivir
trttot=arbidol+10*kaletra+100*tamiflu
data$SPO2[is.na(data$SPO2)]=mean(na.omit(data$SPO2))
x0=cbind(data$Age, data$Sex-1, data$Disease_Status, data$SPO2)
hosp=data$Time2Discharge
####################################
dbfit=function(death, hosp, trt, x){
fitps=glm(trt~x, family="binomial")
ps=fitps$fitted
fit1=glm(death~x, family="binomial", subset=(trt==1))
beta1=fit1$coef
prob1=1/(1+exp(-cbind(1,x)%*%beta1))
mu1death=mean(prob1+trt/ps*(death-prob1))
fit1hosp1=glm((death==0 & hosp<=7)~x, family="binomial", subset=(trt==1))
beta1hosp1=fit1hosp1$coef
prob1hosp1=1/(1+exp(-cbind(1,x)%*%beta1hosp1))
mu1wk1=mean(prob1hosp1+trt/ps*((death==0 & hosp<=7)-prob1hosp1))
fit1hosp2=glm((death==0 & hosp<=14 & hosp>7)~x, family="binomial", subset=(trt==1))
beta1hosp2=fit1hosp2$coef
prob1hosp2=1/(1+exp(-cbind(1,x)%*%beta1hosp2))
mu1wk2=mean(prob1hosp2+trt/ps*((death==0 & hosp<=14 & hosp>7)-prob1hosp2))
mu1wk3=1-mu1wk1-mu1wk2-mu1death
fit0=glm(death~x, family="binomial", subset=(trt==0))
beta0=fit0$coef
prob0=1/(1+exp(-cbind(1,x)%*%beta0))
mu0death=mean(prob0+(1-trt)/(1-ps)*(death-prob0))
fit0hosp1=glm((death==0 & hosp<=7)~x, family="binomial", subset=(trt==0))
beta0hosp1=fit0hosp1$coef
prob0hosp1=1/(1+exp(-cbind(1,x)%*%beta0hosp1))
mu0wk1=mean(prob0hosp1+trt/ps*((death==0 & hosp<=7)-prob0hosp1))
fit0hosp2=glm((death==0 & hosp<=14 & hosp>7)~x, family="binomial", subset=(trt==0))
beta0hosp2=fit0hosp2$coef
prob0hosp2=1/(1+exp(-cbind(1,x)%*%beta0hosp2))
mu0wk2=mean(prob0hosp2+trt/ps*((death==0 & hosp<=14 & hosp>7)-prob0hosp2))
mu0wk3=1-mu0wk1-mu0wk2-mu0death
return(cbind(c(mu1death, mu1wk3, mu1wk2, mu1wk1), c(mu0death, mu0wk3, mu0wk2, mu0wk1)))
}
Compare kaletra
trt=kaletra
x=cbind(x0)
####################################################################
mu1=apply(x0[trt==1,], 2, mean)
mu0=apply(x0[trt==0,], 2, mean)
sd1=apply(x0[trt==1,], 2, sd)
sd0=apply(x0[trt==0,], 2, sd)
pv=rep(0,4)
pv[1]=t.test(x0[trt==1,1], x0[trt==0, 1])$p.value
pv[2]=fisher.test(cbind(c(sum(x0[trt==1, 2]==1), sum(x0[trt==1, 2]==0)), c(sum(x0[trt==0, 2]==1), sum(x0[trt==0, 2]==0))))$p.value
pv[3]=fisher.test(cbind(c(sum(x0[trt==1, 3]==1), sum(x0[trt==1, 3]==0)), c(sum(x0[trt==0, 3]==1), sum(x0[trt==0, 3]==0))))$p.value
pv[4]=t.test(x0[trt==1,4], x0[trt==0, 4])$p.value
balance.table.kaletra=cbind(mu0, sd0, mu1, sd1, pv)
colnames(balance.table.kaletra)=c("mean0", "sd0", "mean1", "sd1", "p-value")
rownames(balance.table.kaletra)=c("age", "sex", "disease", "SPO2")
y=death
p0=mean(y[trt==0])
p1=mean(y[trt==1])
n0=sum(trt==0)
n1=sum(trt==1)
result.kaletra.dth.unadj=cbind( c(p0, p0+c(-1,1)*1.96*sqrt(p0*(1-p0)/n0), NA),
c(p1, p1+c(-1,1)*1.96*sqrt(p1*(1-p1)/n1), NA),
c(p0-p1, p0-p1+c(-1,1)*1.96*sqrt(p0*(1-p0)/n0+p1*(1-p1)/n1), 1-pchisq((p0-p1)^2/(p0*(1-p0)/n0+p1*(1-p1)/n1), 1)))
colnames(result.kaletra.dth.unadj)=c("mortality0", "mortality1", "difference in mortaltiy")
rownames(result.kaletra.dth.unadj)=c("est", "ci-low", "ci-high", "p-value")
y=(death==0 & hosp<=7)
p0=mean(y[trt==0])
p1=mean(y[trt==1])
n0=sum(trt==0)
n1=sum(trt==1)
result.kaletra.wk1.unadj=cbind( c(p0, p0+c(-1,1)*1.96*sqrt(p0*(1-p0)/n0), NA),
c(p1, p1+c(-1,1)*1.96*sqrt(p1*(1-p1)/n1), NA),
c(p0-p1, p0-p1+c(-1,1)*1.96*sqrt(p0*(1-p0)/n0+p1*(1-p1)/n1), 1-pchisq((p0-p1)^2/(p0*(1-p0)/n0+p1*(1-p1)/n1), 1)))
colnames(result.kaletra.wk1.unadj)=c("wk1-arm0", "wk1-arm1", "difference in wk1")
rownames(result.kaletra.wk1.unadj)=c("est", "ci-low", "ci-high", "p-value")
y=(death==0 & hosp<=14 & hosp>7)
p0=mean(y[trt==0])
p1=mean(y[trt==1])
n0=sum(trt==0)
n1=sum(trt==1)
result.kaletra.wk2.unadj=cbind( c(p0, p0+c(-1,1)*1.96*sqrt(p0*(1-p0)/n0), NA),
c(p1, p1+c(-1,1)*1.96*sqrt(p1*(1-p1)/n1), NA),
c(p0-p1, p0-p1+c(-1,1)*1.96*sqrt(p0*(1-p0)/n0+p1*(1-p1)/n1), 1-pchisq((p0-p1)^2/(p0*(1-p0)/n0+p1*(1-p1)/n1), 1)))
colnames(result.kaletra.wk2.unadj)=c("wk2-arm0", "wk2-arm1", "difference in wk2")
rownames(result.kaletra.wk2.unadj)=c("est", "ci-low", "ci-high", "p-value")
y=(death==0 & hosp>14)
p0=mean(y[trt==0])
p1=mean(y[trt==1])
n0=sum(trt==0)
n1=sum(trt==1)
result.kaletra.wk3.unadj=cbind( c(p0, p0+c(-1,1)*1.96*sqrt(p0*(1-p0)/n0), NA),
c(p1, p1+c(-1,1)*1.96*sqrt(p1*(1-p1)/n1), NA),
c(p0-p1, p0-p1+c(-1,1)*1.96*sqrt(p0*(1-p0)/n0+p1*(1-p1)/n1), 1-pchisq((p0-p1)^2/(p0*(1-p0)/n0+p1*(1-p1)/n1), 1)))
colnames(result.kaletra.wk3.unadj)=c("wk3-arm0", "wk3-arm1", "difference in wk3")
rownames(result.kaletra.wk3.unadj)=c("est", "ci-low", "ci-high", "p-value")
Unadjusted summary
kable(round(result.kaletra.dth.unadj, 3))
|
|
mortality0
|
mortality1
|
difference in mortaltiy
|
|
est
|
0.223
|
0.181
|
0.042
|
|
ci-low
|
0.162
|
0.130
|
-0.038
|
|
ci-high
|
0.284
|
0.233
|
0.122
|
|
p-value
|
NA
|
NA
|
0.302
|
kable(round(cbind(result.kaletra.wk1.unadj, result.kaletra.wk2.unadj, result.kaletra.wk3.unadj), 3))
|
|
wk1-arm0
|
wk1-arm1
|
difference in wk1
|
wk2-arm0
|
wk2-arm1
|
difference in wk2
|
wk3-arm0
|
wk3-arm1
|
difference in wk3
|
|
est
|
0.145
|
0.093
|
0.052
|
0.128
|
0.423
|
-0.295
|
0.503
|
0.302
|
0.200
|
|
ci-low
|
0.094
|
0.054
|
-0.012
|
0.079
|
0.357
|
-0.377
|
0.430
|
0.241
|
0.105
|
|
ci-high
|
0.197
|
0.132
|
0.117
|
0.178
|
0.489
|
-0.213
|
0.576
|
0.364
|
0.296
|
|
p-value
|
NA
|
NA
|
0.113
|
NA
|
NA
|
0.000
|
NA
|
NA
|
0.000
|
##Adjusted Analysis (doubly robust regression)
fit=dbfit(death, hosp, trt, x)
mu1=fit[,1]
mu0=fit[,2]
delta=mu0-mu1
n=length(trt)
B=5000
mu1.bt=mu0.bt=matrix(0, B, 4)
for(b in 1:B)
{ set.seed(b)
id=sample(n, n, T)
fit.bt=dbfit(death[id], hosp[id], trt[id], x[id,])
mu1.bt[b,]=fit.bt[,1]
mu0.bt[b,]=fit.bt[,2]
}
delta.bt=mu0.bt-mu1.bt
result.kaletra.dth.adj=cbind(
c(mu0[1], mu0[1]+c(-1, 1)*1.96*sd(mu0.bt[,1]), NA),
c(mu1[1], mu1[1]+c(-1, 1)*1.96*sd(mu1.bt[,1]), NA),
c(delta[1], delta[1]+c(-1,1)*1.96*sd(delta.bt[,1]), 1-pchisq(delta[1]^2/var(delta.bt[,1]),1))
)
colnames(result.kaletra.dth.adj)=c("mortality0", "mortality1", "difference in mortality")
rownames(result.kaletra.dth.adj)=c("est", "ci-low", "ci-high", "p-value")
result.kaletra.wk3.adj=cbind(
c(mu0[2], mu0[2]+c(-1, 1)*1.96*sd(mu0.bt[,2]), NA),
c(mu1[2], mu1[2]+c(-1, 1)*1.96*sd(mu1.bt[,2]), NA),
c(delta[2], delta[2]+c(-1,1)*1.96*sd(delta.bt[,2]), 1-pchisq(delta[2]^2/var(delta.bt[,2]),1))
)
colnames(result.kaletra.wk3.adj)=c("wk1-arm0", "wk1-arm1", "difference in wk1")
rownames(result.kaletra.wk3.adj)=c("est", "ci-low", "ci-high", "p-value")
result.kaletra.wk2.adj=cbind(
c(mu0[3], mu0[3]+c(-1, 1)*1.96*sd(mu0.bt[,3]), NA),
c(mu1[3], mu1[3]+c(-1, 1)*1.96*sd(mu1.bt[,3]), NA),
c(delta[3], delta[3]+c(-1,1)*1.96*sd(delta.bt[,3]), 1-pchisq(delta[3]^2/var(delta.bt[,3]),1))
)
colnames(result.kaletra.wk2.adj)=c("wk2-arm0", "wk2-arm1", "difference in wk2")
rownames(result.kaletra.wk2.adj)=c("est", "ci-low", "ci-high", "p-value")
result.kaletra.wk1.adj=cbind(
c(mu0[4], mu0[4]+c(-1, 1)*1.96*sd(mu0.bt[,4]), NA),
c(mu1[4], mu1[4]+c(-1, 1)*1.96*sd(mu1.bt[,4]), NA),
c(delta[4], delta[4]+c(-1,1)*1.96*sd(delta.bt[,4]), 1-pchisq(delta[4]^2/var(delta.bt[,4]),1))
)
colnames(result.kaletra.wk1.adj)=c("wk3-arm0", "wk3-arm1", "difference in wk3")
rownames(result.kaletra.wk1.adj)=c("est", "ci-low", "ci-high", "p-value")
kable(round(balance.table.kaletra,3))
|
|
mean0
|
sd0
|
mean1
|
sd1
|
p-value
|
|
age
|
58.419
|
15.909
|
60.023
|
13.826
|
0.291
|
|
sex
|
0.419
|
0.495
|
0.470
|
0.500
|
0.360
|
|
disease
|
0.453
|
0.499
|
0.581
|
0.494
|
0.012
|
|
SPO2
|
91.860
|
11.398
|
92.042
|
8.983
|
0.863
|
kable(round(result.kaletra.dth.unadj, 3))
|
|
mortality0
|
mortality1
|
difference in mortaltiy
|
|
est
|
0.223
|
0.181
|
0.042
|
|
ci-low
|
0.162
|
0.130
|
-0.038
|
|
ci-high
|
0.284
|
0.233
|
0.122
|
|
p-value
|
NA
|
NA
|
0.302
|
kable(round(result.kaletra.dth.adj, 3))
|
|
mortality0
|
mortality1
|
difference in mortality
|
|
est
|
0.219
|
0.188
|
0.031
|
|
ci-low
|
0.158
|
0.141
|
-0.039
|
|
ci-high
|
0.279
|
0.235
|
0.102
|
|
p-value
|
NA
|
NA
|
0.384
|
kable(round(cbind(result.kaletra.wk1.unadj, result.kaletra.wk2.unadj, result.kaletra.wk3.unadj), 3))
|
|
wk1-arm0
|
wk1-arm1
|
difference in wk1
|
wk2-arm0
|
wk2-arm1
|
difference in wk2
|
wk3-arm0
|
wk3-arm1
|
difference in wk3
|
|
est
|
0.145
|
0.093
|
0.052
|
0.128
|
0.423
|
-0.295
|
0.503
|
0.302
|
0.200
|
|
ci-low
|
0.094
|
0.054
|
-0.012
|
0.079
|
0.357
|
-0.377
|
0.430
|
0.241
|
0.105
|
|
ci-high
|
0.197
|
0.132
|
0.117
|
0.178
|
0.489
|
-0.213
|
0.576
|
0.364
|
0.296
|
|
p-value
|
NA
|
NA
|
0.113
|
NA
|
NA
|
0.000
|
NA
|
NA
|
0.000
|
kable(round(cbind(result.kaletra.wk1.adj, result.kaletra.wk2.adj, result.kaletra.wk3.adj), 3))
|
|
wk3-arm0
|
wk3-arm1
|
difference in wk3
|
wk2-arm0
|
wk2-arm1
|
difference in wk2
|
wk1-arm0
|
wk1-arm1
|
difference in wk1
|
|
est
|
0.103
|
0.106
|
-0.003
|
0.427
|
0.427
|
0.000
|
0.251
|
0.279
|
-0.028
|
|
ci-low
|
0.063
|
0.063
|
-0.015
|
0.362
|
0.362
|
-0.008
|
0.168
|
0.216
|
-0.101
|
|
ci-high
|
0.143
|
0.149
|
0.009
|
0.492
|
0.492
|
0.008
|
0.334
|
0.343
|
0.044
|
|
p-value
|
NA
|
NA
|
0.617
|
NA
|
NA
|
0.997
|
NA
|
NA
|
0.444
|