---
title: "MSM with IPTW"
author: "(Sukie) Xiuqi Yang"
date: "11/9/2022"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list=ls())
#install packages (if needed)
install.packages("tableone")
install.packages("ipw")
install.packages("sandwich")
install.packages("survey")

#load packages
library(tableone)
library(ipw)
library(sandwich) #for robust variance estimation
library(survey)
library(jtools)
help(ipw)
```

## R Markdown


```{r load data}
expit <- function(x) {1/(1+exp(-x)) }
logit <- function(p) {log(p)-log(1-p)}

#read in data
load(url('https://hbiostat.org/data/repo/rhc.sav'))
#view data
 
View(rhc)

#treatment variables is swang1
#x variables that we will use
#cat1: primary disease category
#age
#sex
#meanbp1: mean blood pressure

#create a data set with just these variables, for simplicity
ARF<-as.numeric(rhc$cat1=='ARF')
CHF<-as.numeric(rhc$cat1=='CHF')
Cirr<-as.numeric(rhc$cat1=='Cirrhosis')
colcan<-as.numeric(rhc$cat1=='Colon Cancer')
Coma<-as.numeric(rhc$cat1=='Coma')
COPD<-as.numeric(rhc$cat1=='COPD')
lungcan<-as.numeric(rhc$cat1=='Lung Cancer')
MOSF<-as.numeric(rhc$cat1=='MOSF w/Malignancy')
sepsis<-as.numeric(rhc$cat1=='MOSF w/Sepsis')
female<-as.numeric(rhc$sex=='Female')
died<-as.integer(rhc$death=='Yes')
age<-rhc$age
treatment<-as.numeric(rhc$swang1=='RHC')
meanbp1<-rhc$meanbp1
#new dataset
mydata<-cbind(ARF,CHF,Cirr,colcan,Coma,lungcan,MOSF,sepsis,
              age,female,meanbp1,treatment,died)
mydata<-data.frame(mydata)

#covariates we will use (shorter list than you would use in practice)
xvars<-c("age","female","meanbp1","ARF","CHF","Cirr","colcan",
         "Coma","lungcan","MOSF","sepsis")


```


```{r table one, echo=FALSE}
#first fit propensity score model to get weights
#look at a table 1
table1<- CreateTableOne(vars=xvars,strata="treatment", data=mydata, test=FALSE)
## include standardized mean difference (SMD)
print(table1,smd=TRUE)

```

```{r mannual, echo=FALSE}

#propensity score model
psmodel <- glm(treatment ~ age + female + meanbp1+ARF+CHF+Cirr+colcan+
         Coma+lungcan+MOSF+sepsis,
               family  = binomial(link ="logit"))

## value of propensity score for each subject
ps <-predict(psmodel, type = "response")

#create weights
weight<-ifelse(treatment==1,1/(ps),1/(1-ps))

#apply weights to data
weighteddata<-svydesign(ids = ~ 1, data =mydata, weights = ~ weight)

#weighted table 1
weightedtable <-svyCreateTableOne(vars = xvars, strata = "treatment", 
                                 data = weighteddata, test = FALSE)
## Show table with SMD
print(weightedtable, smd = TRUE)

## fit the following MSM using IPTW
#get causal risk difference
glm.obj<-glm(died~treatment,weights=weight,family=quasibinomial(link="identity"))
#summary(glm.obj)
betaiptw<-coef(glm.obj)
SE<-sqrt(diag(vcovHC(glm.obj, type="HC0")))

causalrd<-(betaiptw[2])
lcl<-(betaiptw[2]-1.96*SE[2])
ucl<-(betaiptw[2]+1.96*SE[2])
c(lcl,causalrd,ucl)

#get causal relative risk. Weighted GLM
glm.obj<-glm(died~treatment,weights=weight,family=quasibinomial(link=log))
#summary(glm.obj)
betaiptw<-coef(glm.obj)
#to properly account for weighting, use asymptotic (sandwich) variance
SE<-sqrt(diag(vcovHC(glm.obj, type="HC0")))

#get point estimate and CI for relative risk (need to exponentiate)
causalrr<-exp(betaiptw[2])
lcl<-exp(betaiptw[2]-1.96*SE[2])
ucl<-exp(betaiptw[2]+1.96*SE[2])
c(lcl,causalrr,ucl)

#truncate weights at 10

truncweight<-replace(weight,weight>10,10)
#get causal risk difference
glm.obj<-glm(died~treatment,weights=truncweight,family=quasibinomial(link="identity"))
#summary(glm.obj)
betaiptw<-coef(glm.obj)
SE<-sqrt(diag(vcovHC(glm.obj, type="HC0")))

causalrd<-(betaiptw[2])
lcl<-(betaiptw[2]-1.96*SE[2])
ucl<-(betaiptw[2]+1.96*SE[2])
c(lcl,causalrd,ucl)
```

## Written package
```{r ipw_package, echo=FALSE}
# fit propensity score model to get weights, but truncated
weightmodel<-ipwpoint(exposure= treatment, family = "binomial", link ="logit",
                      denominator= ~ age + female + meanbp1+ARF+CHF+Cirr+colcan+
                        Coma+lungcan+MOSF+sepsis, data=mydata,trunc=.01)

#numeric summary of weights
summary(weightmodel$weights.trun)
#plot of weights
ipwplot(weights = weightmodel$weights.trun, logscale = FALSE,
        main = "weights", xlim = c(0, 22))
mydata$wt<-weightmodel$weights.trun
#fit a marginal structural model (risk difference)
msm <- (svyglm(died ~ treatment, design = svydesign(~ 1, weights = ~wt, data =mydata)))
coef(msm)
confint(msm)
summary(msm)
```
## Written package
```{r ipw, echo=FALSE}

mydata <- read.csv("C:/Users/10369/Dropbox/mydata_example.csv")
mydata$parent_outm<-as.factor(mydata$parent_outm)
mydata$type3lag1<-as.factor(mydata$parent_outmlag1)
mydata$p_edu<-as.factor(mydata$p_edu)
mydata$carearr_night<-as.factor(mydata$carearr_night)
mydata$carearr_nightlag1<-as.factor(mydata$carearr_nightlag1)
mydata$lnfinc_pct_5<-as.factor(mydata$lnfinc_pct_5)
mydata$lnfinc_pctlag1_5<-as.factor(mydata$lnfinc_pctlag1_5)

weightmodel1<-ipwtm(exposure= parent_outm,  
                   family = "binomial", link ="logit",
                   numerator= ~sex + ruralorigin + p_edu + bornunderw + BornHosp +familysize + parent_outmlag1, 
                   denominator= ~ sex + ruralorigin + p_edu + bornunderw + BornHosp +familysize + parent_outmlag1 + agemons + lnfinc_pct_5+ lnfinc_pctlag1_5 + carearr_night+ carearr_nightlag1  ,
                   id=pid, timevar = time , 
                   type="all", 
                   data=mydata, maxit=1000)



summary(weightmodel1$ipw.weights)
sd(weightmodel1$ipw.weights)
## also get the standard deviation of it


mydata$wt<-weightmodel1$ipw.weights

msm_noncog <-glm( noncog ~  parent_outm + sex + ruralorigin + p_edu + bornunderw + BornHosp  +familysize , weights = mydata$ipw.weights, data =mydata, gaussian(link = "identity"))
#summary(msm_noncog, digit = 3, robust ="HC3")
summ(msm_noncog, digit = 3, robust ="HC3")
```

