Skip to content

Instantly share code, notes, and snippets.

@ttungl

ttungl/iptw_output

Created June 23, 2022 05:49
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save ttungl/8142a92e11d3b93f5102d03e5cee49ca to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save ttungl/8142a92e11d3b93f5102d03e5cee49ca to your computer and use it in GitHub Desktop.
Download ZIP
IPTW terminal output
Raw
iptw_output
R version 4.2.0 (2022-04-22) -- "Vigorous Calisthenics"
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin20 (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
[Workspace loaded from ~/.RData]
> #install packages (if needed)
> install.packages("tableone")
trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/tableone_0.13.2.tgz'
Content type 'application/x-gzip' length 385087 bytes (376 KB)
==================================================
downloaded 376 KB
The downloaded binary packages are in
/var/folders/5h/pbylctx95sj27h0hdl2b99rw0000gn/T//Rtmp1aKRLU/downloaded_packages
> install.packages("ipw")
trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/ipw_1.0-11.tgz'
Content type 'application/x-gzip' length 196094 bytes (191 KB)
==================================================
downloaded 191 KB
The downloaded binary packages are in
/var/folders/5h/pbylctx95sj27h0hdl2b99rw0000gn/T//Rtmp1aKRLU/downloaded_packages
> install.packages("sandwich")
trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/sandwich_3.0-2.tgz'
Content type 'application/x-gzip' length 1454756 bytes (1.4 MB)
==================================================
downloaded 1.4 MB
The downloaded binary packages are in
/var/folders/5h/pbylctx95sj27h0hdl2b99rw0000gn/T//Rtmp1aKRLU/downloaded_packages
> install.packages("survey")
trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/survey_4.1-1.tgz'
Content type 'application/x-gzip' length 2724306 bytes (2.6 MB)
==================================================
downloaded 2.6 MB
The downloaded binary packages are in
/var/folders/5h/pbylctx95sj27h0hdl2b99rw0000gn/T//Rtmp1aKRLU/downloaded_packages
> library(tableone)
> library(ipw)
> library(sandwich) #for robust variance estimation
> library(survey)
Loading required package: grid
Loading required package: Matrix
Loading required package: survival
Attaching package: ‘survey’
The following object is masked from ‘package:graphics’:
dotchart
> expit <- function(x) {1/(1+exp(-x)) }
> logit <- function(p) {log(p)-log(1-p)}
> load(url("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.sav"))
Error in load(url("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.sav")) :
cannot open the connection to 'http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.sav'
In addition: Warning message:
In load(url("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.sav")) :
URL 'http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.sav': status was 'Couldn't resolve host name'
> load(url("https://biostat.app.vumc.org/wiki/pub/Main/DataSets/rhc.sav"))
> View(rhc)
> #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)
> xvars<-c("age","female","meanbp1","ARF","CHF","Cirr","colcan",
+ "Coma","lungcan","MOSF","sepsis")
> table1<- CreateTableOne(vars=xvars,strata="treatment", data=mydata, test=FALSE)
> ## include standardized mean difference (SMD)
> print(table1,smd=TRUE)
Stratified by treatment
0 1 SMD
n 3551 2184
age (mean (SD)) 61.76 (17.29) 60.75 (15.63) 0.061
female (mean (SD)) 0.46 (0.50) 0.41 (0.49) 0.093
meanbp1 (mean (SD)) 84.87 (38.87) 68.20 (34.24) 0.455
ARF (mean (SD)) 0.45 (0.50) 0.42 (0.49) 0.059
CHF (mean (SD)) 0.07 (0.25) 0.10 (0.29) 0.095
Cirr (mean (SD)) 0.05 (0.22) 0.02 (0.15) 0.145
colcan (mean (SD)) 0.00 (0.04) 0.00 (0.02) 0.038
Coma (mean (SD)) 0.10 (0.29) 0.04 (0.20) 0.207
lungcan (mean (SD)) 0.01 (0.10) 0.00 (0.05) 0.095
MOSF (mean (SD)) 0.07 (0.25) 0.07 (0.26) 0.018
sepsis (mean (SD)) 0.15 (0.36) 0.32 (0.47) 0.415
> #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)
Stratified by treatment
0 1 SMD
n 5732.49 5744.88
age (mean (SD)) 61.36 (17.56) 61.43 (15.33) 0.004
female (mean (SD)) 0.45 (0.50) 0.45 (0.50) 0.001
meanbp1 (mean (SD)) 78.60 (37.58) 79.26 (40.31) 0.017
ARF (mean (SD)) 0.44 (0.50) 0.44 (0.50) 0.010
CHF (mean (SD)) 0.08 (0.27) 0.08 (0.27) 0.005
Cirr (mean (SD)) 0.04 (0.19) 0.04 (0.19) 0.001
colcan (mean (SD)) 0.00 (0.04) 0.00 (0.06) 0.042
Coma (mean (SD)) 0.08 (0.26) 0.07 (0.25) 0.023
lungcan (mean (SD)) 0.01 (0.08) 0.01 (0.09) 0.014
MOSF (mean (SD)) 0.07 (0.26) 0.07 (0.26) 0.004
sepsis (mean (SD)) 0.21 (0.41) 0.22 (0.41) 0.002
> #to get a weighted mean for a single covariate directly:
> mean(weight[treatment==1]*age[treatment==1])/(mean(weight[treatment==1]))
[1] 61.42933
> #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)
treatment treatment treatment
0.02333223 0.05154951 0.07976679
> #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)
treatment treatment treatment
1.036698 1.081764 1.128790
> 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)
treatment treatment treatment
0.02768882 0.05493243 0.08217605
> #first fit propensity score model to get weights
> weightmodel<-ipwpoint(exposure= treatment, family = "binomial", link ="logit",
+ denominator= ~ age + female + meanbp1+ARF+CHF+Cirr+colcan+
+ Coma+lungcan+MOSF+sepsis, data=mydata)
> summary(weightmodel$ipw.weights)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.046 1.405 1.721 2.001 2.280 21.606
> ipwplot(weights = weightmodel$ipw.weights, logscale = FALSE,
+ main = "weights", xlim = c(0, 22))
> mydata$wt<-weightmodel$ipw.weights
> #fit a marginal structural model (risk difference)
> msm <- (svyglm(died ~ treatment, design = svydesign(~ 1, weights = ~wt,
+ data =mydata)))
> coef(msm)
(Intercept) treatment
0.63046375 0.05154951
> confint(msm)
2.5 % 97.5 %
(Intercept) 0.61401097 0.64691653
treatment 0.02332433 0.07977469
> # 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)
> summary(weightmodel$weights.trun)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.081 1.405 1.721 1.972 2.280 6.379
> 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)
(Intercept) treatment
0.63045533 0.05494865
> confint(msm)
2.5 % 97.5 %
(Intercept) 0.61400312 0.64690753
treatment 0.02822366 0.08167363
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment