###############################################################################
#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 tamiflu
trt=tamiflu
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.tamiflu=cbind(mu0, sd0, mu1, sd1, pv)
colnames(balance.table.tamiflu)=c("mean0", "sd0", "mean1", "sd1", "p-value")
rownames(balance.table.tamiflu)=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.tamiflu.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.tamiflu.dth.unadj)=c("mortality0", "mortality1", "difference in mortaltiy")
rownames(result.tamiflu.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.tamiflu.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.tamiflu.wk1.unadj)=c("wk1-arm0", "wk1-arm1", "difference in wk1")
rownames(result.tamiflu.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.tamiflu.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.tamiflu.wk2.unadj)=c("wk2-arm0", "wk2-arm1", "difference in wk2")
rownames(result.tamiflu.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.tamiflu.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.tamiflu.wk3.unadj)=c("wk3-arm0", "wk3-arm1", "difference in wk3")
rownames(result.tamiflu.wk3.unadj)=c("est", "ci-low", "ci-high", "p-value")
Unadjusted summary
kable(round(result.tamiflu.dth.unadj, 3))
|
|
mortality0
|
mortality1
|
difference in mortaltiy
|
|
est
|
0.215
|
0.125
|
0.090
|
|
ci-low
|
0.171
|
0.044
|
-0.002
|
|
ci-high
|
0.259
|
0.206
|
0.183
|
|
p-value
|
NA
|
NA
|
0.056
|
kable(round(cbind(result.tamiflu.wk1.unadj, result.tamiflu.wk2.unadj, result.tamiflu.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.136
|
0.016
|
0.121
|
0.282
|
0.328
|
-0.046
|
0.367
|
0.531
|
-0.165
|
|
ci-low
|
0.099
|
-0.015
|
0.073
|
0.233
|
0.213
|
-0.171
|
0.315
|
0.409
|
-0.297
|
|
ci-high
|
0.173
|
0.046
|
0.169
|
0.330
|
0.443
|
0.079
|
0.419
|
0.654
|
-0.032
|
|
p-value
|
NA
|
NA
|
0.000
|
NA
|
NA
|
0.467
|
NA
|
NA
|
0.015
|
##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.tamiflu.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.tamiflu.dth.adj)=c("mortality0", "mortality1", "difference in mortality")
rownames(result.tamiflu.dth.adj)=c("est", "ci-low", "ci-high", "p-value")
result.tamiflu.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.tamiflu.wk3.adj)=c("wk1-arm0", "wk1-arm1", "difference in wk1")
rownames(result.tamiflu.wk3.adj)=c("est", "ci-low", "ci-high", "p-value")
result.tamiflu.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.tamiflu.wk2.adj)=c("wk2-arm0", "wk2-arm1", "difference in wk2")
rownames(result.tamiflu.wk2.adj)=c("est", "ci-low", "ci-high", "p-value")
result.tamiflu.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.tamiflu.wk1.adj)=c("wk3-arm0", "wk3-arm1", "difference in wk3")
rownames(result.tamiflu.wk1.adj)=c("est", "ci-low", "ci-high", "p-value")
kable(round(balance.table.tamiflu,3))
|
|
mean0
|
sd0
|
mean1
|
sd1
|
p-value
|
|
age
|
59.679
|
15.031
|
57.312
|
13.556
|
0.213
|
|
sex
|
0.436
|
0.497
|
0.500
|
0.504
|
0.410
|
|
disease
|
0.530
|
0.500
|
0.484
|
0.504
|
0.585
|
|
SPO2
|
91.824
|
10.669
|
92.656
|
6.795
|
0.422
|
kable(round(result.tamiflu.dth.unadj, 3))
|
|
mortality0
|
mortality1
|
difference in mortaltiy
|
|
est
|
0.215
|
0.125
|
0.090
|
|
ci-low
|
0.171
|
0.044
|
-0.002
|
|
ci-high
|
0.259
|
0.206
|
0.183
|
|
p-value
|
NA
|
NA
|
0.056
|
kable(round(result.tamiflu.dth.adj, 3))
|
|
mortality0
|
mortality1
|
difference in mortality
|
|
est
|
0.211
|
0.146
|
0.065
|
|
ci-low
|
0.167
|
0.063
|
-0.026
|
|
ci-high
|
0.255
|
0.230
|
0.156
|
|
p-value
|
NA
|
NA
|
0.163
|
kable(round(cbind(result.tamiflu.wk1.unadj, result.tamiflu.wk2.unadj, result.tamiflu.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.136
|
0.016
|
0.121
|
0.282
|
0.328
|
-0.046
|
0.367
|
0.531
|
-0.165
|
|
ci-low
|
0.099
|
-0.015
|
0.073
|
0.233
|
0.213
|
-0.171
|
0.315
|
0.409
|
-0.297
|
|
ci-high
|
0.173
|
0.046
|
0.169
|
0.330
|
0.443
|
0.079
|
0.419
|
0.654
|
-0.032
|
|
p-value
|
NA
|
NA
|
0.000
|
NA
|
NA
|
0.467
|
NA
|
NA
|
0.015
|
kable(round(cbind(result.tamiflu.wk1.adj, result.tamiflu.wk2.adj, result.tamiflu.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.021
|
0.018
|
0.004
|
0.322
|
0.327
|
-0.005
|
0.446
|
0.509
|
-0.063
|
|
ci-low
|
-0.011
|
-0.023
|
-0.030
|
0.210
|
0.210
|
-0.030
|
0.326
|
0.374
|
-0.164
|
|
ci-high
|
0.053
|
0.059
|
0.038
|
0.433
|
0.444
|
0.020
|
0.566
|
0.644
|
0.038
|
|
p-value
|
NA
|
NA
|
0.830
|
NA
|
NA
|
0.677
|
NA
|
NA
|
0.220
|