# redo of week2 Lindner example using optmatch etc # Rogosa R-session, warts and all R version 3.5.3 (2019-03-11) -- "Great Truth" > library(MatchIt) > library(optmatch) Loading required package: survival The optmatch package has an academic license. Enter relaxinfo() for more information. > library(cobalt) > install.packages("PSAgraphics") > library(PSAgraphics) Loading required package: rpart > data(lindner) > head(lindner) lifepres cardbill abcix stent height female diabetic acutemi ejecfrac ves1proc 1 0.0 14301 1 0 163 1 1 0 56 1 2 11.6 3563 1 0 168 0 0 0 56 1 3 11.6 4694 1 0 188 0 0 0 50 1 4 11.6 7366 1 0 175 0 1 0 50 1 5 11.6 8247 1 0 168 1 0 0 55 1 6 11.6 8319 1 0 178 0 0 0 50 1 > covs = subset(lindner, select = -c(lifepres, cardbill, abcix )) > head(covs) stent height female diabetic acutemi ejecfrac ves1proc 1 0 163 1 1 0 56 1 2 0 168 0 0 0 56 1 3 0 188 0 0 0 50 1 4 0 175 0 1 0 50 1 5 0 168 1 0 0 55 1 6 0 178 0 0 0 50 1 #start the optmatch process > pfit = glm(f.build("abcix", covs), data = lindner, family = "binomial") > lindner$p.score = pfit$fitted.values #get the propensity score > boxplot(lindner$p.score ~ lindner$abcix) # propensity score > fm2 = fullmatch(abcix ~ p.score, data = lindner) > bal.tab(fm2, formula = f.build("abcix", covs), data = lindner) Call fullmatch(x = abcix ~ p.score, data = lindner) Balance Measures Type Diff.Adj stent Binary -0.0185 height Contin. -0.0099 female Binary 0.0347 diabetic Binary -0.0136 acutemi Binary 0.0124 ejecfrac Contin. -0.0570 ves1proc Contin. -0.0020 Sample sizes Control Treated All 298 698 Matched 298 698 > # gives you spectacular balance > #optmatch balance much better on stent female diabetic > love.plot(bal.tab(fm2, formula = f.build("abcix", covs), data = lindner), threshold = .1) > # outcome analysis log(cardbill) > # look at subclasses > summary(fm2) Structure of matched sets: 5+:1 4:1 3:1 2:1 1:1 1:2 1:3 1:4 1:5+ 38 14 34 40 113 14 4 2 2 Effective Sample Size: 337.8 (equivalent number of matched pairs). > stratumStructure(fm2) 17:1 16:1 15:1 13:1 12:1 11:1 10:1 9:1 8:1 7:1 6:1 5:1 4:1 3:1 2:1 1:1 1:2 1:3 1:4 1:5 1:6 3 2 1 1 2 1 1 2 4 3 8 10 14 34 40 113 14 4 2 1 1 > # looks like 261 subclasses compared to 267 from MatchIt(full) > # here we grab a factor giving us the subclass info > # I call the augmented lindner dataset "Lmatched" > > Lmatched = cbind(lindner, matches = fm2) > head(Lmatched) lifepres cardbill abcix stent height female diabetic acutemi ejecfrac ves1proc p.score matches 1 0.0 14301 1 0 163 1 1 0 56 1 0.4079170 1.1 2 11.6 3563 1 0 168 0 0 0 56 1 0.5784602 1.112 3 11.6 4694 1 0 188 0 0 0 50 1 0.5244469 1.223 4 11.6 7366 1 0 175 0 1 0 50 1 0.4727311 1.334 5 11.6 8247 1 0 168 1 0 0 55 1 0.4930466 1.445 6 11.6 8319 1 0 178 0 0 0 50 1 0.5625524 1.556 > str(Lmatched) 'data.frame': 996 obs. of 12 variables: $ lifepres: num 0 11.6 11.6 11.6 11.6 11.6 11.6 11.6 11.6 11.6 ... $ cardbill: int 14301 3563 4694 7366 8247 8319 8410 8517 8763 8823 ... $ abcix : int 1 1 1 1 1 1 1 1 1 1 ... $ stent : int 0 0 0 0 0 0 0 0 0 0 ... $ height : int 163 168 188 175 168 178 185 173 152 180 ... $ female : int 1 0 0 0 1 0 0 1 1 0 ... $ diabetic: int 1 0 0 1 0 0 0 0 0 0 ... $ acutemi : int 0 0 0 0 0 0 0 0 0 0 ... $ ejecfrac: int 56 56 50 50 55 50 58 30 60 60 ... $ ves1proc: int 1 1 1 1 1 1 1 1 1 1 ... $ p.score : num 0.408 0.578 0.524 0.473 0.493 ... $ matches : Factor w/ 261 levels "1.1","1.10","1.101",..: 1 8 84 142 182 213 230 236 241 70 ... ..- attr(*, "exceedances")= Named num 0.565 .. ..- attr(*, "names")= chr "1" ..- attr(*, "call")= language fullmatch(x = abcix ~ p.score, data = lindner) ..- attr(*, "contrast.group")= logi TRUE TRUE TRUE TRUE TRUE TRUE ... ..- attr(*, "subproblem")= Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ... .. ..- attr(*, "names")= chr "1" "2" "3" "4" ... ..- attr(*, "min.controls")= num 0 ..- attr(*, "max.controls")= num Inf ..- attr(*, "omit.fraction")= num NA ..- attr(*, "hashed.distance")= chr "3fb101a69d79c551c6755deb7ec872dc" > length(table(Lmatched$matches)) [1] 261 > table(Lmatched$matches) # lots of small subclasses, causes me lmer problems 1.1 1.10 1.101 1.104 1.106 1.107 1.110 1.112 1.113 1.119 1.12 1.120 1.122 1.124 1.125 1.126 1.129 1.13 1.131 1.134 1.135 2 2 2 2 4 2 2 8 3 2 4 11 2 6 2 2 3 2 2 2 3 1.136 1.14 1.143 1.144 1.147 1.148 1.149 1.15 1.150 1.152 1.153 1.154 1.155 1.157 1.158 1.16 1.161 1.162 1.163 1.164 1.166 3 4 2 2 3 4 2 2 2 5 3 3 4 4 7 2 3 4 6 2 6 1.167 1.168 1.171 1.172 1.173 1.174 1.175 1.176 1.178 1.179 1.18 1.180 1.181 1.182 1.183 1.184 1.185 1.188 1.189 1.19 1.190 2 2 2 3 2 5 2 5 2 3 4 4 4 5 4 7 2 3 2 9 2 1.191 1.192 1.193 1.194 1.196 1.197 1.2 1.200 1.201 1.204 1.206 1.209 1.21 1.212 1.213 1.215 1.216 1.219 1.22 1.221 1.223 4 4 2 2 2 2 2 2 2 4 3 3 2 5 5 2 2 5 5 3 6 1.229 1.231 1.233 1.234 1.235 1.238 1.24 1.244 1.246 1.25 1.251 1.253 1.254 1.257 1.26 1.261 1.263 1.264 1.268 1.27 1.272 3 2 2 2 2 2 6 3 2 3 2 2 2 2 2 2 3 3 3 6 2 1.275 1.279 1.28 1.285 1.287 1.29 1.292 1.298 1.299 1.3 1.30 1.300 1.301 1.302 1.303 1.304 1.305 1.306 1.307 1.308 1.31 3 4 3 4 2 2 4 4 18 3 3 5 2 2 5 2 6 3 9 3 3 1.310 1.311 1.312 1.314 1.315 1.317 1.318 1.319 1.32 1.323 1.327 1.328 1.329 1.33 1.330 1.334 1.335 1.336 1.338 1.34 1.340 3 4 3 17 6 7 8 4 4 6 6 2 2 2 3 3 2 3 6 3 7 1.343 1.35 1.351 1.353 1.354 1.357 1.36 1.361 1.363 1.364 1.368 1.37 1.377 1.379 1.38 1.385 1.39 1.397 1.398 1.4 1.40 2 4 7 9 2 2 2 7 3 2 2 2 2 2 3 4 5 2 2 2 4 1.401 1.408 1.41 1.410 1.412 1.415 1.418 1.422 1.423 1.43 1.431 1.437 1.44 1.445 1.457 1.46 1.468 1.469 1.470 1.471 1.473 3 2 4 4 2 3 3 2 4 5 2 3 2 5 2 3 3 18 8 17 13 1.474 1.477 1.478 1.48 1.480 1.483 1.485 1.486 1.490 1.495 1.5 1.50 1.503 1.512 1.517 1.519 1.52 1.524 1.53 1.535 1.54 10 5 4 2 13 18 16 12 2 10 2 2 14 4 9 3 3 7 2 3 4 1.545 1.55 1.556 1.557 1.560 1.562 1.57 1.579 1.58 1.59 1.601 1.61 1.623 1.63 1.634 1.64 1.65 1.656 1.66 1.666 1.667 2 4 2 3 4 7 3 2 2 3 2 3 2 3 3 2 2 2 4 2 5 1.669 1.67 1.673 1.676 1.677 1.679 1.68 1.686 1.687 1.688 1.690 1.692 1.695 1.697 1.7 1.70 1.73 1.74 1.75 1.76 1.79 2 4 2 7 2 2 2 3 4 4 2 3 2 2 2 2 2 2 2 5 3 1.8 1.81 1.83 1.84 1.88 1.89 1.9 1.90 1.94 2 3 4 2 2 2 2 4 2 > library(lme4) # lmer doesn't like me today > optmatch_lmer = lmer(log(cardbill) ~ abcix +(abcix|matches), data = Lmatched) boundary (singular) fit: see ?isSingular > ?isSingular starting httpd help server ... done > isSingular(optmatch_lmer, tol = 1e-03) [1] TRUE > isSingular(optmatch_lmer, tol = 1e-02) [1] TRUE > isSingular(optmatch_lmer, tol = 1e-06) [1] FALSE # appease lmer with a control statement > optmatch_lmer = lmer(log(cardbill) ~ abcix +(abcix|matches), control = lmerControl(check.conv.singular = .makeCC(action = "message", tol = 1e-6)) ,data = Lmatched) > summary(optmatch_lmer) Linear mixed model fit by REML ['lmerMod'] Formula: log(cardbill) ~ abcix + (abcix | matches) Data: Lmatched Control: lmerControl(check.conv.singular = .makeCC(action = "message", tol = 1e-06)) REML criterion at convergence: 1278.1 Scaled residuals: Min 1Q Median 3Q Max -3.2461 -0.6354 -0.2558 0.4198 4.3152 Random effects: Groups Name Variance Std.Dev. Corr matches (Intercept) 0.08324 0.2885 abcix 0.07030 0.2651 -1.00 Residual 0.18820 0.4338 Number of obs: 996, groups: matches, 261 Fixed effects: Estimate Std. Error t value (Intercept) 9.40325 0.03132 300.201 abcix 0.17531 0.03471 5.051 Correlation of Fixed Effects: (Intr) abcix -0.879 > # about the same as Matchit lmer > confint(optmatch_lmer, oldNames = FALSE) Computing profile confidence intervals ... 2.5 % 97.5 % sd_(Intercept)|matches 0.1914737 0.3702391 cor_abcix.(Intercept)|matches -1.0000000 1.0000000 sd_abcix|matches 0.1595840 0.3665216 sigma 0.4093783 0.4568357 (Intercept) 9.3418121 9.4649082 abcix 0.1070127 0.2441981 There were 12 warnings (use warnings() to see them) > # lmer is appeased, about same result; lmer from Matchit may well now have same issue