###############################################################################
#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 arbidol
trt=arbidol
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.arbidol=cbind(mu0, sd0, mu1, sd1, pv)
colnames(balance.table.arbidol)=c("mean0", "sd0", "mean1", "sd1", "p-value")
rownames(balance.table.arbidol)=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.arbidol.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.arbidol.dth.unadj)=c("mortality0", "mortality1", "difference in mortaltiy")
rownames(result.arbidol.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.arbidol.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.arbidol.wk1.unadj)=c("wk1-arm0", "wk1-arm1", "difference in wk1")
rownames(result.arbidol.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.arbidol.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.arbidol.wk2.unadj)=c("wk2-arm0", "wk2-arm1", "difference in wk2")
rownames(result.arbidol.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.arbidol.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.arbidol.wk3.unadj)=c("wk3-arm0", "wk3-arm1", "difference in wk3")
rownames(result.arbidol.wk3.unadj)=c("est", "ci-low", "ci-high", "p-value")
Unadjusted summary
kable(round(result.arbidol.dth.unadj, 3))
|
|
mortality0
|
mortality1
|
difference in mortaltiy
|
|
est
|
0.310
|
0.091
|
0.218
|
|
ci-low
|
0.245
|
0.051
|
0.142
|
|
ci-high
|
0.374
|
0.132
|
0.294
|
|
p-value
|
NA
|
NA
|
0.000
|
kable(round(cbind(result.arbidol.wk1.unadj, result.arbidol.wk2.unadj, result.arbidol.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.168
|
0.066
|
0.102
|
0.406
|
0.173
|
0.234
|
0.117
|
0.670
|
-0.553
|
|
ci-low
|
0.115
|
0.031
|
0.039
|
0.338
|
0.120
|
0.147
|
0.072
|
0.604
|
-0.633
|
|
ci-high
|
0.220
|
0.101
|
0.164
|
0.475
|
0.225
|
0.320
|
0.162
|
0.736
|
-0.474
|
|
p-value
|
NA
|
NA
|
0.001
|
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.arbidol.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.arbidol.dth.adj)=c("mortality0", "mortality1", "difference in mortality")
rownames(result.arbidol.dth.adj)=c("est", "ci-low", "ci-high", "p-value")
result.arbidol.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.arbidol.wk3.adj)=c("wk1-arm0", "wk1-arm1", "difference in wk1")
rownames(result.arbidol.wk3.adj)=c("est", "ci-low", "ci-high", "p-value")
result.arbidol.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.arbidol.wk2.adj)=c("wk2-arm0", "wk2-arm1", "difference in wk2")
rownames(result.arbidol.wk2.adj)=c("est", "ci-low", "ci-high", "p-value")
result.arbidol.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.arbidol.wk1.adj)=c("wk3-arm0", "wk3-arm1", "difference in wk3")
rownames(result.arbidol.wk1.adj)=c("est", "ci-low", "ci-high", "p-value")
kable(round(balance.table.arbidol,3))
|
|
mean0
|
sd0
|
mean1
|
sd1
|
p-value
|
|
age
|
59.325
|
14.975
|
59.264
|
14.684
|
0.968
|
|
sex
|
0.447
|
0.498
|
0.447
|
0.498
|
1.000
|
|
disease
|
0.538
|
0.500
|
0.508
|
0.501
|
0.614
|
|
SPO2
|
90.370
|
13.177
|
93.548
|
5.234
|
0.002
|
kable(round(result.arbidol.dth.unadj, 3))
|
|
mortality0
|
mortality1
|
difference in mortaltiy
|
|
est
|
0.310
|
0.091
|
0.218
|
|
ci-low
|
0.245
|
0.051
|
0.142
|
|
ci-high
|
0.374
|
0.132
|
0.294
|
|
p-value
|
NA
|
NA
|
0.000
|
kable(round(result.arbidol.dth.adj, 3))
|
|
mortality0
|
mortality1
|
difference in mortality
|
|
est
|
0.282
|
0.101
|
0.181
|
|
ci-low
|
0.224
|
0.057
|
0.112
|
|
ci-high
|
0.340
|
0.145
|
0.251
|
|
p-value
|
NA
|
NA
|
0.000
|
kable(round(cbind(result.arbidol.wk1.unadj, result.arbidol.wk2.unadj, result.arbidol.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.168
|
0.066
|
0.102
|
0.406
|
0.173
|
0.234
|
0.117
|
0.670
|
-0.553
|
|
ci-low
|
0.115
|
0.031
|
0.039
|
0.338
|
0.120
|
0.147
|
0.072
|
0.604
|
-0.633
|
|
ci-high
|
0.220
|
0.101
|
0.164
|
0.475
|
0.225
|
0.320
|
0.162
|
0.736
|
-0.474
|
|
p-value
|
NA
|
NA
|
0.001
|
NA
|
NA
|
0.000
|
NA
|
NA
|
0.000
|
kable(round(cbind(result.arbidol.wk1.adj, result.arbidol.wk2.adj, result.arbidol.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.073
|
0.064
|
0.009
|
0.174
|
0.173
|
0.001
|
0.471
|
0.662
|
-0.191
|
|
ci-low
|
0.039
|
0.031
|
-0.002
|
0.123
|
0.121
|
-0.009
|
0.391
|
0.593
|
-0.262
|
|
ci-high
|
0.107
|
0.098
|
0.020
|
0.225
|
0.225
|
0.011
|
0.551
|
0.731
|
-0.120
|
|
p-value
|
NA
|
NA
|
0.111
|
NA
|
NA
|
0.863
|
NA
|
NA
|
0.000
|
Analysis for two drug combination (only use PS method )
trttot=arbidol+10*kaletra
table(trttot)
## trttot
## 0 1 10 11
## 68 111 129 86
### 0: no medicine
### 1: arbidol only
### 10: kaletra only
### 11: arbidol+kaletra
psfit=function(death, trttot, x0)
{
fit0=glm((trttot==0)~x0, family="binomial")
fit1=glm((trttot==1)~x0, family="binomial")
fit10=glm((trttot==10)~x0, family="binomial")
fit11=glm((trttot==11)~x0, family="binomial")
ps0=fit0$fitted.values
ps1=fit1$fitted.values
ps10=fit10$fitted.values
ps11=fit11$fitted.values
mu0=mean((trttot==0)/ps0*death)/mean((trttot==0)/ps0)
mu1=mean((trttot==1)/ps1*death)/mean((trttot==1)/ps1)
mu10=mean((trttot==10)/ps10*death)/mean((trttot==10)/ps10)
mu11=mean((trttot==11)/ps11*death)/mean((trttot==11)/ps11)
return(c(mu0, mu1, mu10, mu11))
}
muest=psfit(death, trttot, x0)
n=length(trttot)
B=1000
mubt=matrix(NA, B, 4)
for(b in 1:B)
{ set.seed(b)
id=sample(n, n, T)
mubt[b,]=psfit(death[id], trttot[id], x0[id,])
}
result=cbind(muest, muest-1.96*apply(mubt,2,sd), muest+1.96*apply(mubt,2, sd))
colnames(result)=c("est", "ci-low", "ci-high")
rownames(result)=c("none", "arbidol only ", "kaletra only", "arbidol+kaletra")
kable(round(result,3))
|
|
est
|
ci-low
|
ci-high
|
|
none
|
0.353
|
0.252
|
0.453
|
|
arbidol only
|
0.096
|
0.033
|
0.159
|
|
kaletra only
|
0.231
|
0.164
|
0.299
|
|
arbidol+kaletra
|
0.091
|
0.036
|
0.146
|
May not have enough sample size for three drug combinations
trttot=arbidol+10*kaletra+100*tamiflu
table(trttot)
## trttot
## 0 1 10 11 100 101 110 111
## 66 92 99 73 2 19 30 13
fit0=glm((trttot==0)~x0, family="binomial")
fit1=glm((trttot==1)~x0, family="binomial")
fit10=glm((trttot==10)~x0, family="binomial")
fit11=glm((trttot==11)~x0, family="binomial")
fit110=glm((trttot==110)~x0, family="binomial")
fit111=glm((trttot==111)~x0, family="binomial")
ps0=fit0$fitted.values
ps1=fit1$fitted.values
ps10=fit10$fitted.values
ps11=fit11$fitted.values
ps110=fit110$fitted.values
ps111=fit111$fitted.values
mu0=mean((trttot==0)/ps0*death)
mu1=mean((trttot==1)/ps1*death)
mu10=mean((trttot==10)/ps10*death)
mu11=mean((trttot==11)/ps11*death)
mu110=mean((trttot==110)/ps110*death)
mu111=mean((trttot==111)/ps111*death)
### 0: no medicine
### 1: arbidol only
### 10: kaletra only
### 11: arbidol+kaletra
### 110: kaletra+tamiflu
### 111: arbidol+kaletra+tamiflu