# install.packages('optmatch') install.packages('lme4')
# install.packages('MatchIt') install.packages('plyr') library(knitr)
library(MatchIt)
## Loading required package: MASS
library(optmatch)
## Loading required package: digest You're loading optmatch, by B. Hansen and
## M. Fredrickson. The optmatch package makes essential use of D. P.
## Bertsekas and P. Tseng's RELAX-IV algorithm and code, as well as
## Bertsekas' AUCTION algorithm and code. Using the software to 'satisfy in
## any part commercial delivery requirements to government or industry'
## requires a special agreement with Dr. Bertsekas. For more information,
## enter relaxinfo() at the command line.
library(lme4)
## Loading required package: lattice Loading required package: Matrix
library(plyr)
rm(list=ls())
setwd("C:/Users/Yan/Dropbox/Projects_Friends/DIF_Propensity_Score/analysis/Yan/try_out")
timss <- foreign::read.spss("can_imp1_book1_Dec17.sav", to.data.frame = TRUE)
## Warning: can_imp1_book1_Dec17.sav: Unrecognized record type 7, subtype 18
## encountered in system file Warning: can_imp1_book1_Dec17.sav: Unrecognized
## record type 7, subtype 24 encountered in system file
str(timss[1:40])
## 'data.frame': 822 obs. of 40 variables:
## $ IDCNTRY : Factor w/ 3 levels " 9132"," 9133",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ IDSCHOOL: Factor w/ 448 levels " 189"," 190",..: 1 2 3 3 4 5 6 7 7 8 ...
## $ IDCLASS : Factor w/ 545 levels " 18901"," 19004",..: 1 2 3 3 4 5 6 7 8 9 ...
## $ IDSTUD : Factor w/ 822 levels " 1890110"," 1900413",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ BSDAGE : num 14 14 14 14 14 14 15 14 13 15 ...
## $ ITSEX : Factor w/ 2 levels "BOY","GIRL": 1 2 2 2 2 1 1 2 2 2 ...
## $ ITLANG : Factor w/ 2 levels "French","Engish": 2 2 2 2 2 2 2 2 2 2 ...
## $ BS4GBOOK: num 3 3 2 2 4 2 3 3 4 4 ...
## $ BS4MHCAL: num 3 3 3 3 3 3 3 2 2 3 ...
## $ BS4GUSIN: num 2 4 3 3 0 3 3 2 2 2 ...
## $ BSDGEDUP: num 4 4 2 3 4 4 4 3 4 4 ...
## $ BSDGCAVL: num 4 3 4 3 3 3 3 4 4 4 ...
## $ BSDMTMH : num 1 1 1 1 2 2 1 1 1 1 ...
## $ BSDMPATM: num 2 2 0 0 2 2 0 2 2 1 ...
## $ BSDMSVM : num 2 2 2 2 2 2 2 2 2 2 ...
## $ BSDMSCM : num 2 2 0 1 2 1 2 2 2 2 ...
## $ BSDGPBSS: num 2 0 1 0 1 1 1 2 2 2 ...
## $ BSMMAT01: num 581 529 620 559 613 462 522 625 550 585 ...
## $ BC4GTENR: num 719 606 512 512 373 327 266 203 203 186 ...
## $ BC4GCOMU: num 5 5 5 5 5 5 3 5 5 1 ...
## $ BC4GSBED: num 1 0 0 0 0 0 0 0 0 0 ...
## $ BC4GSBEA: num 1 2 3 3 2 2 2 3 3 3 ...
## $ BC4GNALA: num 2 3 2 2 3 3 3 3 3 1 ...
## $ BC4GCHTS: num 4 3 3 3 4 3 2 3 3 3 ...
## $ BC4GCHES: num 2 2 3 3 3 1 2 2 2 2 ...
## $ BC4GCHPS: num 3 3 3 3 3 4 3 4 4 3 ...
## $ BC4GCHPI: num 4 1 2 2 2 2 2 4 4 0 ...
## $ BC4MSOEM: Factor w/ 2 levels "NO","YES": 1 2 1 1 2 2 2 2 2 1 ...
## $ BC4MSORM: Factor w/ 2 levels "NO","YES": 1 2 1 1 2 2 1 2 2 2 ...
## $ BCDGAS : num 2 2 1 1 2 1 1 1 1 1 ...
## $ BCDSRMI : num 2 2 2 2 2 1 1 2 2 2 ...
## $ BCDGPPSC: num 2 1 1 1 2 2 1 2 2 1 ...
## $ BCDGHW : num 25 29 25 25 28 25 37 27 27 28 ...
## $ Y1 : num 1 0 0 1 1 1 1 1 1 1 ...
## $ Y2 : num 1 1 1 1 1 0 0 1 1 1 ...
## $ Y3 : num 0 1 1 1 1 1 0 1 0 0 ...
## $ Y4 : num 0 1 0 0 0 0 0 0 0 0 ...
## $ Y5 : num 1 0 1 1 1 0 0 0 1 0 ...
## $ Y6 : num 1 0 0 1 1 1 1 1 1 1 ...
## $ Y7 : num 1 1 1 1 1 0 1 1 1 1 ...
summary(timss)
## IDCNTRY IDSCHOOL IDCLASS IDSTUD BSDAGE
## 9132:250 257: 5 318604: 4 1890110: 1 Min. :11
## 9133:271 3252: 5 25702: 3 1900413: 1 1st Qu.:14
## 9135:301 3253: 5 320301: 3 1910107: 1 Median :14
## 3270: 5 321605: 3 1910125: 1 Mean :14
## 3274: 5 321708: 3 1930214: 1 3rd Qu.:14
## 249: 4 323308: 3 1950103: 1 Max. :17
## (Other) :793 (Other) :803 (Other) :816
## ITSEX ITLANG BS4GBOOK BS4MHCAL BS4GUSIN
## BOY :378 French:281 Min. :0.00 Min. :0.00 Min. :0.00
## GIRL:444 Engish:541 1st Qu.:1.00 1st Qu.:2.00 1st Qu.:1.00
## Median :2.00 Median :3.00 Median :2.00
## Mean :2.16 Mean :2.25 Mean :2.06
## 3rd Qu.:3.00 3rd Qu.:3.00 3rd Qu.:3.00
## Max. :4.00 Max. :3.00 Max. :4.00
##
## BSDGEDUP BSDGCAVL BSDMTMH BSDMPATM
## Min. :0.00 Min. :0.00 Min. :0.00 Min. :0.00
## 1st Qu.:3.00 1st Qu.:3.00 1st Qu.:1.00 1st Qu.:0.00
## Median :4.00 Median :4.00 Median :1.00 Median :1.00
## Mean :3.23 Mean :3.65 Mean :1.13 Mean :1.14
## 3rd Qu.:4.00 3rd Qu.:4.00 3rd Qu.:2.00 3rd Qu.:2.00
## Max. :4.00 Max. :4.00 Max. :2.00 Max. :2.00
##
## BSDMSVM BSDMSCM BSDGPBSS BSMMAT01 BC4GTENR
## Min. :0.00 Min. :0.00 Min. :0.00 Min. :278 Min. : 20
## 1st Qu.:2.00 1st Qu.:1.00 1st Qu.:1.00 1st Qu.:467 1st Qu.: 432
## Median :2.00 Median :2.00 Median :2.00 Median :516 Median : 724
## Mean :1.78 Mean :1.39 Mean :1.41 Mean :516 Mean : 859
## 3rd Qu.:2.00 3rd Qu.:2.00 3rd Qu.:2.00 3rd Qu.:567 3rd Qu.:1224
## Max. :2.00 Max. :2.00 Max. :2.00 Max. :734 Max. :2400
##
## BC4GCOMU BC4GSBED BC4GSBEA BC4GNALA
## Min. :0.0 Min. :0.000 Min. :0.00 Min. :0.00
## 1st Qu.:2.0 1st Qu.:0.000 1st Qu.:0.00 1st Qu.:1.00
## Median :4.0 Median :1.000 Median :1.00 Median :3.00
## Mean :3.2 Mean :0.989 Mean :1.28 Mean :2.04
## 3rd Qu.:5.0 3rd Qu.:2.000 3rd Qu.:2.00 3rd Qu.:3.00
## Max. :5.0 Max. :3.000 Max. :3.00 Max. :3.00
##
## BC4GCHTS BC4GCHES BC4GCHPS BC4GCHPI BC4MSOEM
## Min. :0.00 Min. :0 Min. :0.00 Min. :0.00 NO :435
## 1st Qu.:2.00 1st Qu.:2 1st Qu.:2.00 1st Qu.:1.00 YES:387
## Median :3.00 Median :2 Median :2.00 Median :2.00
## Mean :2.79 Mean :2 Mean :2.43 Mean :1.78
## 3rd Qu.:3.00 3rd Qu.:2 3rd Qu.:3.00 3rd Qu.:2.00
## Max. :4.00 Max. :3 Max. :4.00 Max. :4.00
##
## BC4MSORM BCDGAS BCDSRMI BCDGPPSC BCDGHW
## NO :146 Min. :0.000 Min. :0.00 Min. :0.00 Min. :15.0
## YES:676 1st Qu.:1.000 1st Qu.:1.00 1st Qu.:1.00 1st Qu.:25.0
## Median :1.000 Median :1.00 Median :1.00 Median :25.0
## Mean :0.971 Mean :1.46 Mean :1.19 Mean :25.8
## 3rd Qu.:1.000 3rd Qu.:2.00 3rd Qu.:2.00 3rd Qu.:26.0
## Max. :2.000 Max. :2.00 Max. :2.00 Max. :45.0
##
## Y1 Y2 Y3 Y4
## Min. :0.000 Min. :0.000 Min. :0.00 Min. :0.000
## 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:0.00 1st Qu.:0.000
## Median :1.000 Median :1.000 Median :1.00 Median :0.000
## Mean :0.781 Mean :0.816 Mean :0.63 Mean :0.337
## 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:1.00 3rd Qu.:1.000
## Max. :1.000 Max. :1.000 Max. :1.00 Max. :1.000
##
## Y5 Y6 Y7 Y8
## Min. :0.00 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:0.00 1st Qu.:0.250 1st Qu.:0.000 1st Qu.:0.000
## Median :0.00 Median :1.000 Median :1.000 Median :0.000
## Mean :0.42 Mean :0.749 Mean :0.631 Mean :0.378
## 3rd Qu.:1.00 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:1.000
## Max. :1.00 Max. :1.000 Max. :1.000 Max. :1.000
##
## Y9 Y10 Y11 Y12
## Min. :0.000 Min. :0.00 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:0.00 1st Qu.:0.000 1st Qu.:0.000
## Median :1.000 Median :0.00 Median :0.000 Median :0.000
## Mean :0.579 Mean :0.38 Mean :0.748 Mean :0.566
## 3rd Qu.:1.000 3rd Qu.:0.00 3rd Qu.:2.000 3rd Qu.:2.000
## Max. :1.000 Max. :2.00 Max. :2.000 Max. :2.000
##
## Y13 Y14 Y15 Y16
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:0.000
## Median :0.000 Median :1.000 Median :1.000 Median :1.000
## Mean :0.461 Mean :0.787 Mean :0.848 Mean :0.593
## 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:1.000
## Max. :1.000 Max. :1.000 Max. :1.000 Max. :1.000
##
## Y17 Y18 Y19 Y20
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.000
## Median :0.000 Median :1.000 Median :1.000 Median :0.000
## Mean :0.466 Mean :0.562 Mean :0.569 Mean :0.478
## 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:1.000
## Max. :1.000 Max. :1.000 Max. :1.000 Max. :1.000
##
## Y21 Y22 Y23 Y24
## Min. :0.000 Min. :0.00 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:0.00 1st Qu.:0.000 1st Qu.:0.000
## Median :0.000 Median :0.00 Median :0.000 Median :1.000
## Mean :0.359 Mean :0.17 Mean :0.377 Mean :0.596
## 3rd Qu.:1.000 3rd Qu.:0.00 3rd Qu.:1.000 3rd Qu.:1.000
## Max. :1.000 Max. :1.00 Max. :1.000 Max. :1.000
##
## Y25 Y26 Y27 Y28
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:1.000 1st Qu.:1.000
## Median :0.000 Median :1.000 Median :1.000 Median :1.000
## Mean :0.456 Mean :0.711 Mean :0.882 Mean :0.854
## 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:1.000
## Max. :1.000 Max. :1.000 Max. :1.000 Max. :1.000
##
## Y29 tot1 tot2 tot3
## Min. :0.00 Min. : 3.0 Min. : 3.0 Min. : 3.0
## 1st Qu.:0.00 1st Qu.:11.0 1st Qu.:11.0 1st Qu.:12.0
## Median :1.00 Median :16.0 Median :16.0 Median :16.0
## Mean :1.05 Mean :16.5 Mean :16.4 Mean :16.6
## 3rd Qu.:2.00 3rd Qu.:21.0 3rd Qu.:21.0 3rd Qu.:21.0
## Max. :2.00 Max. :32.0 Max. :32.0 Max. :32.0
##
## tot4 tot5 tot6 tot7
## Min. : 3.0 Min. : 3.0 Min. : 3.0 Min. : 3.0
## 1st Qu.:12.0 1st Qu.:12.0 1st Qu.:11.0 1st Qu.:12.0
## Median :17.0 Median :16.5 Median :16.0 Median :16.0
## Mean :16.9 Mean :16.8 Mean :16.5 Mean :16.6
## 3rd Qu.:22.0 3rd Qu.:21.8 3rd Qu.:21.0 3rd Qu.:21.0
## Max. :32.0 Max. :32.0 Max. :32.0 Max. :32.0
##
## tot8 tot9 tot10 tot11
## Min. : 3.0 Min. : 3.0 Min. : 3.0 Min. : 3.0
## 1st Qu.:12.0 1st Qu.:12.0 1st Qu.:12.0 1st Qu.:12.0
## Median :16.5 Median :16.0 Median :17.0 Median :16.0
## Mean :16.9 Mean :16.7 Mean :16.9 Mean :16.5
## 3rd Qu.:22.0 3rd Qu.:21.0 3rd Qu.:22.0 3rd Qu.:21.0
## Max. :32.0 Max. :32.0 Max. :31.0 Max. :31.0
##
## tot12 tot13 tot14 tot15
## Min. : 3.0 Min. : 3.0 Min. : 3.0 Min. : 2.0
## 1st Qu.:12.0 1st Qu.:12.0 1st Qu.:11.0 1st Qu.:11.0
## Median :17.0 Median :16.0 Median :16.0 Median :16.0
## Mean :16.7 Mean :16.8 Mean :16.4 Mean :16.4
## 3rd Qu.:21.0 3rd Qu.:21.0 3rd Qu.:21.0 3rd Qu.:21.0
## Max. :31.0 Max. :32.0 Max. :32.0 Max. :32.0
##
## tot16 tot17 tot18 tot19
## Min. : 3.0 Min. : 3.0 Min. : 3.0 Min. : 3.0
## 1st Qu.:12.0 1st Qu.:12.0 1st Qu.:12.0 1st Qu.:12.0
## Median :16.0 Median :16.5 Median :16.0 Median :16.0
## Mean :16.6 Mean :16.8 Mean :16.7 Mean :16.7
## 3rd Qu.:21.0 3rd Qu.:22.0 3rd Qu.:21.0 3rd Qu.:21.0
## Max. :32.0 Max. :32.0 Max. :32.0 Max. :32.0
##
## tot20 tot21 tot22 tot23
## Min. : 3.0 Min. : 3.0 Min. : 3.0 Min. : 3.0
## 1st Qu.:12.0 1st Qu.:12.0 1st Qu.:12.0 1st Qu.:12.0
## Median :16.5 Median :17.0 Median :17.0 Median :17.0
## Mean :16.8 Mean :16.9 Mean :17.1 Mean :16.9
## 3rd Qu.:21.0 3rd Qu.:22.0 3rd Qu.:22.0 3rd Qu.:22.0
## Max. :32.0 Max. :32.0 Max. :32.0 Max. :32.0
##
## tot24 tot25 tot26 tot27
## Min. : 2.0 Min. : 3.0 Min. : 2.0 Min. : 2.0
## 1st Qu.:12.0 1st Qu.:12.0 1st Qu.:12.0 1st Qu.:11.0
## Median :16.0 Median :16.0 Median :16.0 Median :16.0
## Mean :16.6 Mean :16.8 Mean :16.5 Mean :16.4
## 3rd Qu.:21.0 3rd Qu.:21.0 3rd Qu.:21.0 3rd Qu.:21.0
## Max. :32.0 Max. :32.0 Max. :32.0 Max. :32.0
##
## tot28 tot29
## Min. : 3.0 Min. : 3.0
## 1st Qu.:11.0 1st Qu.:12.0
## Median :16.0 Median :16.0
## Mean :16.4 Mean :16.2
## 3rd Qu.:21.0 3rd Qu.:20.0
## Max. :32.0 Max. :31.0
##
######################################### stratification propensity score matching 1-level model example stage one:
######################################### Obtain prpns
psmod <- glm(ITLANG ~ ITSEX + BS4GBOOK + BS4MHCAL + BSDGEDUP + BSDGCAVL + BSDMTMH +
BSDMPATM + BSDMSVM + BSDMSCM + BSDGPBSS, data = timss, family = binomial)
summary(psmod)
##
## Call:
## glm(formula = ITLANG ~ ITSEX + BS4GBOOK + BS4MHCAL + BSDGEDUP +
## BSDGCAVL + BSDMTMH + BSDMPATM + BSDMSVM + BSDMSCM + BSDGPBSS,
## family = binomial, data = timss)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.538 -1.088 0.606 0.880 1.838
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0879 0.6916 1.57 0.1157
## ITSEXGIRL 0.1559 0.1612 0.97 0.3334
## BS4GBOOK 0.4990 0.0702 7.11 1.2e-12 ***
## BS4MHCAL -0.4901 0.0979 -5.01 5.6e-07 ***
## BSDGEDUP -0.2010 0.0911 -2.21 0.0273 *
## BSDGCAVL -0.0034 0.1356 -0.03 0.9800
## BSDMTMH 0.5272 0.1242 4.24 2.2e-05 ***
## BSDMPATM -0.2803 0.1075 -2.61 0.0091 **
## BSDMSVM 0.0724 0.1693 0.43 0.6691
## BSDMSCM 0.1107 0.1234 0.90 0.3698
## BSDGPBSS -0.1938 0.1197 -1.62 0.1053
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1055.87 on 821 degrees of freedom
## Residual deviance: 942.07 on 811 degrees of freedom
## AIC: 964.1
##
## Number of Fisher Scoring iterations: 4
prpns <- fitted(psmod, response = T)
densityplot(prpns)