propensity score script from "Propensity score based data analysis" by Susanne Stampf

nonrandom.R

## Propensity score based data analysis by Susanne Stampf - http://cran.r-project.org/web/packages/nonrandom/vignettes/nonrandom.pdf


#install.packages("nonrandom")
library("nonrandom")

## PS ESTIMATION
require(glm2)
data(pride)
pride.effect  <-  relative.effect(data = pride, sel = c(2:14), resp = 15, treat = "PCR_RSV") # family = binomial
pride.effect

data(stu1)
stu1.effect <- relative.effect(data = stu1, formula = pst~therapie+tgr+age)

stu1.effect

stu1.ps <- pscore(data = stu1, formula = therapie~tgr+age)

pride.ps <- pscore(data = pride, formula = PCR_RSV~SEX+RSVINF+REGION+ AGE+ELTATOP+EINZ+EXT, name.pscore = "ps")

## Figure 1
plot.pscore(pride.ps, main = "PRI.De study", with.legend = TRUE, par.1 = list(lty=1, lwd=2), par.0 = list(lty=3, lwd=2), xlab = "", ylim = c(0,4.5))



## STRATIFICATION

## (1) stratification (factorization)
stu1.str4 <- ps.makestrata(object = stu1.ps)
stu1.str4
## (2) stratification (specification of the number of strata)
pride.str.b3 <- ps.makestrata(object = pride.ps, breaks = 3)
pride.str.b3

## (3) stratification (definition of specific stratum limits)
pride.str5 <- ps.makestrata(object = pride.ps, breaks = quantile(pride.ps$pscore, seq(0,1,0.2)), name.stratum.index = "stratum")





## MATCHING

pride.m1 <- ps.match(object = pride.ps, ratio = 1, x = 0.2, caliper = "logit", matched.by = "ps", setseed = 38902)
summary(pride.m1)
stu1.m2 <- ps.match(object = stu1.ps, ratio = 2, caliper = 0.5, givenTmatchingC = FALSE, setseed = 39062)
stu1.m2




## BALANCE CHECK


## Figure 2 (left)
stu1.dplot1 <- dist.plot(object = stu1.m2, sel = c("tgr"), compare =TRUE, label.match = c("original data", "matched sample"))



## Figure 2 (right)
stu1.dplot2 <- dist.plot(object = stu1.m2, sel = c("age"), plot.type = 2, compare = TRUE, legend.tile = "Therapy")



## Figure 3 (left)
pride.dplot1 <- dist.plot(object = pride.str5, sel = c("AGE"), label.stratum = c("Str.", "Original"))



## Figure 3 (right)

pride.dplot2 <- dist.plot(object = pride.m1, sel = c("AGE"), plot.type = 2, compare = TRUE, legend.title = "RSV infection")



## Return values of 'stu1.dplot1' presenting in Figure 2 (left)
names(stu1.dplot1)
stu1.dplot1$var.cat
stu1.dplot1$frequency



## Return values of 'stu1.dplot2' presenting in Figure 2 (right)
stu1.dplot2$var.cat
stu1.dplot2$x.s.cat
stu1.dplot2$y.s.cat




## Return values of 'pride.dplot1' presenting in Figure 3 (left)
pride.dplot1$var.noncat
pride.dplot1$mean
pride.dplot1$var.noncat
pride.dplot2$breaks.noncat
pride.dplot2$x.s.noncat
pride.dplot2$y.s.noncat



## PRI.De, stratification
## Balance check using statistical tests
pride.str5.bal <- ps.balance(object = pride.str5, sel = c(2:8), cat.levels = 4, method = "classical", alpha = 5)
pride.str5.bal




## STU1, Matching
## Balance check using standardized differences
stu1.m2.bal <- ps.balance(object  = stu1.m2, sel = c("tgr","age"), method  = "stand.diff", alpha   =  20)
stu1.m2.bal

## Illustration of standardized differences, Figure 4 (right)

stu1.distplot <- plot.stdf(x = stu1.m2.bal, sel = c("tgr", "age"), main  = "STU1: Matching by PS", pch.p = c(20,10))




## PROPENSITY SCORE BASED TREATMENT EFFECTS

## STU1: Effect estimation of therapy ('ther') on quality
## of life ('pst') based on PS stratification
stu1.est <- ps.estimate(object  = stu1.str4, resp = "pst", weights = "opt", regr = pst ~therapie + tgr + age)

## PS stratification in STU1 data
## Summary for treatment effect estimation
summary(stu1.est)

## PS stratification in PRI.De
## Effect estimation of exposure to RSV infection on the
## severity of LRTI
pride.est <- ps.estimate(object  = pride.str5, family  = "binomial", resp = "SEVERE", treat = "PCR_RSV", adj = c("REGION", "ETHNO", "AGE"), weights = "rr")
summary(pride.est)

## STU1, matching by PS, continuous outcome 'pst'
## Effect estimation of therapy on quality of life ('pst')
stu1.est.m2 <- ps.estimate(object = stu1.m2, resp = "pst")
stu1.est.m2

## PRI.De, matching on PS, binary outcome 'SEVERE'
## Effect estimation of exposure to RSV infection on the
## severity of LRTI
pride.estimate.m <- ps.estimate(object = pride.m1, resp = "SEVERE", family = "binomial")
pride.estimate.m



# Leituras

* [The Importance of Covariate Selection in Controlling for Selection Bias in Observational Studies](http://cds.web.unc.edu/files/2012/03/Steiner_et-al_2010.pdf)