###############################################################################
#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")
death=data$Death
arbidol=data$Arbidol
kaletra=data$Kaletra
tamiflu=data$Tamiflu
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")round(result.arbidol.dth.unadj, 3)## mortality0 mortality1 difference in mortaltiy
## est 0.297 0.095 0.201
## ci-low 0.225 0.055 0.119
## ci-high 0.369 0.136 0.284
## p-value NA NA 0.000round(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
## est 0.194 NA NA 0.368 NA
## ci-low 0.131 NA NA 0.292 NA
## ci-high 0.256 NA NA 0.444 NA
## p-value NA NA NA NA NA
## difference in wk2 wk3-arm0 wk3-arm1 difference in wk3
## est NA 0.142 NA NA
## ci-low NA 0.087 NA NA
## ci-high NA 0.197 NA NA
## p-value NA NA NA NAAdjusted 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.btresult.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 | 58.665 | 15.109 | 59.588 | 14.844 | 0.566 |
| sex | 0.465 | 0.500 | 0.447 | 0.498 | 0.748 |
| disease | 0.535 | 0.500 | 0.508 | 0.501 | 0.668 |
| SPO2 | 89.993 | 14.047 | 93.477 | 5.377 | 0.004 |
kable(round(result.arbidol.dth.unadj, 3))| mortality0 | mortality1 | difference in mortaltiy | |
|---|---|---|---|
| est | 0.297 | 0.095 | 0.201 |
| ci-low | 0.225 | 0.055 | 0.119 |
| ci-high | 0.369 | 0.136 | 0.284 |
| p-value | NA | NA | 0.000 |
kable(round(result.arbidol.dth.adj, 3))| mortality0 | mortality1 | difference in mortality | |
|---|---|---|---|
| est | 0.264 | 0.098 | 0.166 |
| ci-low | 0.199 | 0.056 | 0.089 |
| ci-high | 0.329 | 0.140 | 0.243 |
| 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.194 | NA | NA | 0.368 | NA | NA | 0.142 | NA | NA |
| ci-low | 0.131 | NA | NA | 0.292 | NA | NA | 0.087 | NA | NA |
| ci-high | 0.256 | NA | NA | 0.444 | NA | NA | 0.197 | NA | NA |
| p-value | NA | NA | NA | NA | NA | NA | NA | NA | NA |
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 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| ci-low | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| ci-high | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| p-value | NA | NA | NA | NA | NA | NA | NA | NA | NA |
Compare Kaletra (NEEDS to be updated as above)
COmpare Tamiflu
Analysis for two drug combination (only use PS method
trttot=arbidol+10*kaletra
table(trttot)## trttot
## 0 1 10 11
## 60 115 95 84### 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.280 | 0.183 | 0.377 |
| arbidol only | 0.097 | 0.038 | 0.157 |
| kaletra only | 0.240 | 0.158 | 0.322 |
| arbidol+kaletra | 0.085 | 0.032 | 0.138 |
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
## 58 96 65 71 2 19 30 13fit0=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