library(nonrandom)
## Warning: package 'nonrandom' was built under R version 3.3.3
## Loading required package: lme4
## Warning: package 'lme4' was built under R version 3.3.3
## Loading required package: Matrix
ca.dat <- read.csv(file.choose(), header=T)
head(ca.dat)
## klinik idnr tmass therapie alter tgr age ewb pst
## 1 3 782 6 1 42 1 1 63.46154 81.25
## 2 3 975 9 1 49 1 1 90.38462 93.75
## 3 4 448 10 1 39 1 1 73.07692 93.75
## 4 6 835 6 1 53 1 1 75.00000 81.25
## 5 6 866 10 1 52 1 1 34.61538 56.25
## 6 6 953 10 1 51 1 1 63.46154 93.75
### Estimate effect of confounders
ca.dat.effect <- relative.effect(pst~therapie+tmass+alter+tgr+age+ewb, treat= therapie, ca.dat)
ca.dat.effect
##
## Treatment: therapie
## Outcome: pst
## Covariates: tmass alter tgr age ewb
##
##
## Unadjusted treatment effect: 1.5894
##
## Adjusted and relative effects:
##
## adj. treatment effect rel. effect
## alter 0.2010374 87.351652
## age 0.7880392 50.420198
## tmass 1.7729920 11.548500
## tgr 1.7004732 6.985956
## ewb 1.6513514 3.895436
### Determine propensity score
ca.dat.ps <- pscore (therapie ~ tmass+alter+tgr+age+ewb, ca.dat)
plot(ca.dat.ps, with.legend=T, xlim=c(0,1))
## Loading required package: colorspace
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'colorspace'
ca.dat.ps
##
## Call: glm(formula = formula, family = family, data = data, na.action = na.action)
##
## Coefficients:
## (Intercept) tmass alter tgr age
## 5.307503 -0.028233 -0.056772 -0.240628 -0.048594
## ewb
## -0.002856
##
## Degrees of Freedom: 645 Total (i.e. Null); 640 Residual
## Null Deviance: 738.4
## Residual Deviance: 681.4 AIC: 693.4
### Stratification by PS
ca.dat.str.b3 <- ps.makestrata(ca.dat.ps, breaks=quantile(ca.dat.ps$pscore))
ca.dat.str.b3
##
## Stratified by: pscore
##
## Strata information:
##
## Strata bounds n n (proportion)
## 1 [0.367,0.643] 162 0.251
## 2 (0.643,0.763] 161 0.249
## 3 (0.763,0.844] 161 0.249
## 4 (0.844,0.968] 162 0.251
### Look at data
names(ca.dat.str.b3)
## [1] "data" "pscore" "name.pscore"
## [4] "formula.pscore" "model.pscore" "treat"
## [7] "name.treat" "stratum.index" "name.stratum.index"
## [10] "intervals" "stratified.by"
head(ca.dat.str.b3$data)
## klinik idnr tmass therapie alter tgr age ewb pst pscore
## 1 3 782 6 1 42 1 1 63.46154 81.25 0.9074784
## 2 3 975 9 1 49 1 1 90.38462 93.75 0.8486744
## 3 4 448 10 1 39 1 1 73.07692 93.75 0.9099607
## 4 6 835 6 1 53 1 1 75.00000 81.25 0.8355911
## 5 6 866 10 1 52 1 1 34.61538 56.25 0.8435579
## 6 6 953 10 1 51 1 1 63.46154 93.75 0.8401499
## stratum.index
## 1 4
## 2 4
## 3 4
## 4 3
## 5 3
## 6 3
### Generate matched pairs
ca.dat.m2 <- ps.match(object= ca.dat.ps, ratio=1:1,
x=0.2, caliper= "logit", givenTmatchingC=F)
## Argument 'givenTmatchingC'=FALSE: Treated elements were matched to each untreated element.
ca.dat.m2
##
## Matched by: pscore
##
## Matching parameter:
##
## Caliper size: 0.147
## Ratio: 1.000
## Who is treated?: 1.000
##
## Matching information:
##
## Untreated to treated?: FALSE
## Best match?: TRUE
##
## Matching data:
##
## Number of treated obs.: 479
## Number of matched treated obs.: 167
## Number of untreated obs.: 167
## Number of matched untreated obs.: 167
## Number of total matched obs.: 334
## Number of not matched obs.: 312
## Number of matching sets: 167
## Number of incomplete matching sets: 0
head(ca.dat.m2$data.matched)
## klinik idnr tmass therapie alter tgr age ewb pst pscore
## 8 7 592 7 1 53 1 1 48.07692 75.00 0.8421655
## 14 10 107 8 1 42 1 1 69.23077 81.25 0.9011695
## 15 11 514 8 0 47 1 1 75.00000 100.00 0.8710136
## 18 15 188 10 0 43 1 1 92.30769 100.00 0.8840292
## 19 16 505 10 0 36 1 1 75.00000 93.75 0.9225837
## 21 18 235 6 1 52 1 1 51.92308 68.75 0.8517575
## match.index
## 8 44
## 14 15
## 15 1
## 18 2
## 19 3
## 21 4
head(ca.dat.m2$data.matched[order(ca.dat.m2$data.matched$match.index),])
## klinik idnr tmass therapie alter tgr age ewb pst pscore
## 15 11 514 8 0 47 1 1 75.00000 100.00 0.8710136
## 150 17 362 14 1 42 2 1 30.76923 56.25 0.8710291
## 18 15 188 10 0 43 1 1 92.30769 100.00 0.8840292
## 49 38 218 10 1 44 1 1 71.15385 87.50 0.8844011
## 19 16 505 10 0 36 1 1 75.00000 93.75 0.9225837
## 255 43 298 20 1 27 2 1 71.15385 87.50 0.9225105
## match.index
## 15 1
## 150 1
## 18 2
## 49 2
## 19 3
## 255 3
### Balance check (stratified)
dplot1 <- dist.plot(ca.dat.str.b3,sel= "alter", compare= T, label.match = c("Original", "Matched"), col= c("dark blue", "grey"), main= "Balance Check")
### Balance check (matched pairs)
dplot1 <- dist.plot(ca.dat.m2,sel= "alter", compare= T, label.match= c("Original", "Matched"),col= c("dark blue", "grey"), main= "Balance Check")
### Perform statistical test in case of matched data
bal.table <- ps.balance(object= ca.dat.m2, sel= c("tgr","alter", "tmass", "age"), alpha=5)
bal.table
##
## Summary of balance check:
##
## Before: no bal (0) Before: bal (1)
## After: no bal (0) 0 0
## After: bal (1) 4 0
##
##
## Covariates not completely tested: ---
##
##
## Detailed balance check (overall):
##
## tgr alter tmass age
## Before 0 0 0 0
## After 1 1 1 1
##
##
## Detailed balance check:
## [p-values from tests (significance level: 0.05)]
##
## tgr alter tmass age
## Before 0.024 0 0.006 0
## ------ ----- ----- ----- -----
## After 0.555 0.73 0.689 0.908
##
## --------- ---- ---- ---- ----
## Test chi^2 t t chi^2
### stratified analysis using regression
mod.str <- lm(pst~therapie+factor(stratum.index), data= ca.dat.str.b3$data)
summary(mod.str)
##
## Call:
## lm(formula = pst ~ therapie + factor(stratum.index), data = ca.dat.str.b3$data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -62.501 -7.963 0.816 10.787 23.287
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 76.7128 1.3388 57.300 <2e-16 ***
## therapie 0.8292 1.3062 0.635 0.5258
## factor(stratum.index)2 3.7092 1.5591 2.379 0.0176 *
## factor(stratum.index)3 2.8920 1.5833 1.827 0.0682 .
## factor(stratum.index)4 4.1294 1.6043 2.574 0.0103 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.98 on 641 degrees of freedom
## Multiple R-squared: 0.01513, Adjusted R-squared: 0.008981
## F-statistic: 2.461 on 4 and 641 DF, p-value: 0.0442
mod.est <- ps.estimate(object= ca.dat.str.b3$data, resp="pst", treat= "therapie", stratum.index= "stratum.index", regr= pst~therapie+stratum.index)
mod.est
##
## Effect estimation for treatment/exposure on outcome
##
## Treatment/exposure: therapie
## Outcome: pst
## Effect measure: difference ('effect')
##
##
## Table of effect estimates:
##
## effect SE[effect] [95%-CI[effect]]
## ------ ---------- ----------------
## Crude 1.589 1.261 [-0.883,4.061]
## Stratification
## Unadjusted 0.23 1.4056 [-2.525,2.985]
## Adjusted [,]
## Regression 0.829 1.3062 [-1.731,3.389]
##
### check unadjusted effect
delta <- c()
for(i in unique(ca.dat.str.b3$data$stratum.index)){
delta <- c(delta, mean(ca.dat.str.b3$data$pst[ca.dat.str.b3$data$therapie == 1 &
ca.dat.str.b3$data$stratum.index == i])-
mean(ca.dat.str.b3$data$pst[ca.dat.str.b3$data$therapie == 0 &
ca.dat.str.b3$data$stratum.index == i]))
}
mean(delta)
## [1] 0.2321721
### Adjusted analysis using regression
mod.adj <- lm(pst~therapie+pscore, ca.dat.ps$data)
summary(mod.adj)
##
## Call:
## lm(formula = pst ~ therapie + pscore, data = ca.dat.ps$data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -61.034 -8.842 1.229 11.032 23.483
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 71.5922 3.2096 22.306 <2e-16 ***
## therapie 0.6886 1.3115 0.525 0.5997
## pscore 10.6618 4.4513 2.395 0.0169 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.98 on 643 degrees of freedom
## Multiple R-squared: 0.01128, Adjusted R-squared: 0.008207
## F-statistic: 3.669 on 2 and 643 DF, p-value: 0.02604
### TABLE OF EFFECT ESTIMATES###
est <- ps.estimate(object= ca.dat.str.b3, resp= "pst", regr= pst~therapie+pscore)
est
##
## Effect estimation for treatment/exposure on outcome
##
## Treatment/exposure: therapie
## Outcome: pst
## Effect measure: difference ('effect')
##
##
## Table of effect estimates:
##
## effect SE[effect] [95%-CI[effect]]
## ------ ---------- ----------------
## Crude 1.589 1.261 [-0.883,4.061]
## Stratification
## Unadjusted 0.23 1.4056 [-2.525,2.985]
## Adjusted [,]
## Regression 0.689 1.3115 [-1.881,3.259]
##