Yan — Sep 29, 2013, 8:46 AM
#install.packages("MatchIt")
#install.packages("optmatch")
#install.packages("lme4")
#library(MatchIt)
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(knitr)
library(markdown)
#rm(timss)
setwd("C:/Users/Yan/Dropbox/Projects_Friends/DIF_Propensity_Score/analysis/Yan/try_out")
#timss<-foreign::read.spss("qat_imp_book1_tot.sav", to.data.frame=TRUE)
timss<-foreign::read.spss("merge_all_book1.sav", to.data.frame=TRUE)
Warning: merge_all_book1.sav: Unrecognized record type 7, subtype 18
encountered in system file Warning: merge_all_book1.sav: Unrecognized
record type 7, subtype 24 encountered in system file
timss.can<-foreign::read.spss("can_imp_book1_tot.sav", to.data.frame=TRUE)
Warning: can_imp_book1_tot.sav: Unrecognized record type 7, subtype 18
encountered in system file Warning: can_imp_book1_tot.sav: Unrecognized
record type 7, subtype 24 encountered in system file
str(timss)
'data.frame': 2354 obs. of 84 variables:
$ IDCNTRY : Factor w/ 5 levels " 48"," 90",..: 1 1 1 1 1 1 1 1 1 1 ...
$ IDSCH : Factor w/ 948 levels " 1"," 2",..: 1 1 1 1 2 2 2 2 2 3 ...
$ IDSCHOOL: Factor w/ 667 levels " 1"," 2",..: 1 1 1 1 2 2 2 2 2 3 ...
$ IDSTUD : Factor w/ 2332 levels " 10108"," 10126",..: 1 2 3 6 9 10 14 15 17 18 ...
$ BC4GSBED: num 1 1 1 1 1 1 1 1 1 1 ...
$ BC4GSBEA: num 3 3 3 3 0 0 0 0 0 2 ...
$ BC4GCHTS: num 1 1 1 1 0 0 0 0 0 1 ...
$ BC4GCHES: num 1 1 1 1 1 1 1 1 1 1 ...
$ BC4GCHPS: num 2 2 2 2 2 2 2 2 2 1 ...
$ BCDGAS : num 2 2 2 2 1 1 1 1 1 1 ...
$ BCDSRMI : num 1 1 1 1 1 1 1 1 1 1 ...
$ BCDGPPSC: num 1 1 1 1 0 0 0 0 0 1 ...
$ ITLANG : Factor w/ 2 levels " 0"," 1": 1 1 1 1 1 1 1 1 1 1 ...
$ ITSEX : Factor w/ 2 levels " 0"," 1": 1 1 1 1 1 1 1 1 1 1 ...
$ BS4GBOOK: num 1 1 1 2 2 0 3 2 2 0 ...
$ BS4MHCAL: num 2 3 3 3 2 2 2 2 2 2 ...
$ BSDGEDUP: num 2 2 2 2 3 0 0 0 2 1 ...
$ BSDGCAVL: num 0 1 4 1 2 0 1 0 0 0 ...
$ BSDMPATM: num 0 0 0 0 0 1 0 0 0 2 ...
$ BSDMSVM : num 0 0 0 0 0 0 0 0 0 0 ...
$ BSDMSCM : num 0 1 0 1 0 2 0 0 0 0 ...
$ BSDGPBSS: num 0 2 1 0 0 0 1 0 0 2 ...
$ BC4GTENR: num 570 570 570 570 480 480 480 480 480 408 ...
$ BCDGHW : num 25 25 25 25 28 28 28 28 28 27 ...
$ BSDAGE : num 14 13 14 14 13 16 14 13 14 14 ...
$ BSMMAT05: num 361 382 352 390 432 265 494 586 401 429 ...
$ m1 : num 1 1 0 0 1 0 1 1 0 1 ...
$ m2 : num 1 0 0 0 0 0 1 1 0 0 ...
$ m3 : num 0 1 0 0 0 1 0 1 1 0 ...
$ m4 : num 0 0 0 1 0 0 0 1 0 0 ...
$ m5 : num 0 0 0 1 0 0 1 1 0 0 ...
$ m6 : num 1 1 0 0 1 1 0 1 1 1 ...
$ m7 : num 0 0 0 1 0 0 0 1 1 0 ...
$ m8 : num 1 0 0 1 0 0 1 1 1 0 ...
$ m9 : num 0 0 0 0 0 0 1 1 0 0 ...
$ m10 : num 0 0 0 0 0 0 0 0 0 0 ...
$ m11 : num 0 0 0 0 0 0 0 0 0 0 ...
$ m12 : num 0 0 0 0 0 0 1 2 0 0 ...
$ m13 : num 0 0 0 0 0 0 1 1 0 0 ...
$ m14 : num 0 1 0 1 0 0 1 1 0 0 ...
$ m15 : num 1 1 0 0 0 0 1 1 0 1 ...
$ m16 : num 0 1 0 0 1 0 1 0 0 0 ...
$ m17 : num 0 0 0 0 1 0 0 1 1 1 ...
$ m18 : num 0 0 0 0 0 0 1 0 1 0 ...
$ m19 : num 0 0 1 0 0 0 1 1 0 1 ...
$ m20 : num 1 0 0 1 0 0 0 1 0 0 ...
$ m21 : num 1 0 0 0 0 0 1 1 0 0 ...
$ m22 : num 0 0 0 0 0 0 1 1 0 0 ...
$ m23 : num 0 0 0 0 0 0 0 0 0 0 ...
$ m24 : num 0 1 1 0 1 0 0 1 0 1 ...
$ m25 : num 1 0 0 0 1 0 1 1 0 0 ...
$ m26 : num 0 1 1 1 1 1 1 1 0 1 ...
$ m27 : num 1 1 1 1 1 0 1 1 0 0 ...
$ m28 : num 0 1 1 0 1 0 0 1 0 1 ...
$ m29 : num 0 2 0 0 0 0 0 2 0 0 ...
$ tot1 : num 8 11 5 8 8 3 16 25 6 7 ...
$ tot2 : num 8 12 5 8 9 3 16 25 6 8 ...
$ tot3 : num 9 11 5 8 9 2 17 25 5 8 ...
$ tot4 : num 9 12 5 7 9 3 17 25 6 8 ...
$ tot5 : num 9 12 5 7 9 3 16 25 6 8 ...
$ tot6 : num 8 11 5 8 8 2 17 25 5 7 ...
$ tot7 : num 9 12 5 7 9 3 17 25 5 8 ...
$ tot8 : num 8 12 5 7 9 3 16 25 5 8 ...
$ tot9 : num 9 12 5 8 9 3 16 25 6 8 ...
$ tot10 : num 9 12 5 8 9 3 17 26 6 8 ...
$ tot11 : num 9 12 5 8 9 3 17 26 6 8 ...
$ tot12 : num 9 12 5 8 9 3 16 24 6 8 ...
$ tot13 : num 9 12 5 8 9 3 16 25 6 8 ...
$ tot14 : num 9 11 5 7 9 3 16 25 6 8 ...
$ tot15 : num 8 11 5 8 9 3 16 25 6 7 ...
$ tot16 : num 9 11 5 8 8 3 16 26 6 8 ...
$ tot17 : num 9 12 5 8 8 3 17 25 5 7 ...
$ tot18 : num 9 12 5 8 9 3 16 26 5 8 ...
$ tot19 : num 9 12 4 8 9 3 16 25 6 7 ...
$ tot20 : num 8 12 5 7 9 3 17 25 6 8 ...
$ tot21 : num 8 12 5 8 9 3 16 25 6 8 ...
$ tot22 : num 9 12 5 8 9 3 16 25 6 8 ...
$ tot23 : num 9 12 5 8 9 3 17 26 6 8 ...
$ tot24 : num 9 11 4 8 8 3 17 25 6 7 ...
$ tot25 : num 8 12 5 8 8 3 16 25 6 8 ...
$ tot26 : num 9 11 4 7 8 2 16 25 6 7 ...
$ tot27 : num 8 11 4 7 8 3 16 25 6 8 ...
$ tot28 : num 9 11 4 8 8 3 17 25 6 7 ...
$ tot29 : num 9 10 5 8 9 3 17 24 6 8 ...
- attr(*, "codepage")= int 65001
head(timss)
IDCNTRY IDSCH IDSCHOOL IDSTUD BC4GSBED BC4GSBEA BC4GCHTS BC4GCHES
1 48 1 1 10108 1 3 1 1
2 48 1 1 10126 1 3 1 1
3 48 1 1 10314 1 3 1 1
4 48 1 1 10520 1 3 1 1
5 48 2 2 20208 1 0 0 1
6 48 2 2 20226 1 0 0 1
BC4GCHPS BCDGAS BCDSRMI BCDGPPSC ITLANG ITSEX BS4GBOOK BS4MHCAL
1 2 2 1 1 0 0 1 2
2 2 2 1 1 0 0 1 3
3 2 2 1 1 0 0 1 3
4 2 2 1 1 0 0 2 3
5 2 1 1 0 0 0 2 2
6 2 1 1 0 0 0 0 2
BSDGEDUP BSDGCAVL BSDMPATM BSDMSVM BSDMSCM BSDGPBSS BC4GTENR BCDGHW
1 2 0 0 0 0 0 570 25
2 2 1 0 0 1 2 570 25
3 2 4 0 0 0 1 570 25
4 2 1 0 0 1 0 570 25
5 3 2 0 0 0 0 480 28
6 0 0 1 0 2 0 480 28
BSDAGE BSMMAT05 m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14 m15 m16
1 14 361 1 1 0 0 0 1 0 1 0 0 0 0 0 0 1 0
2 13 382 1 0 1 0 0 1 0 0 0 0 0 0 0 1 1 1
3 14 352 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
4 14 390 0 0 0 1 1 0 1 1 0 0 0 0 0 1 0 0
5 13 432 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1
6 16 265 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0
m17 m18 m19 m20 m21 m22 m23 m24 m25 m26 m27 m28 m29 tot1 tot2 tot3 tot4
1 0 0 0 1 1 0 0 0 1 0 1 0 0 8 8 9 9
2 0 0 0 0 0 0 0 1 0 1 1 1 2 11 12 11 12
3 0 0 1 0 0 0 0 1 0 1 1 1 0 5 5 5 5
4 0 0 0 1 0 0 0 0 0 1 1 0 0 8 8 8 7
5 1 0 0 0 0 0 0 1 1 1 1 1 0 8 9 9 9
6 0 0 0 0 0 0 0 0 0 1 0 0 0 3 3 2 3
tot5 tot6 tot7 tot8 tot9 tot10 tot11 tot12 tot13 tot14 tot15 tot16 tot17
1 9 8 9 8 9 9 9 9 9 9 8 9 9
2 12 11 12 12 12 12 12 12 12 11 11 11 12
3 5 5 5 5 5 5 5 5 5 5 5 5 5
4 7 8 7 7 8 8 8 8 8 7 8 8 8
5 9 8 9 9 9 9 9 9 9 9 9 8 8
6 3 2 3 3 3 3 3 3 3 3 3 3 3
tot18 tot19 tot20 tot21 tot22 tot23 tot24 tot25 tot26 tot27 tot28 tot29
1 9 9 8 8 9 9 9 8 9 8 9 9
2 12 12 12 12 12 12 11 12 11 11 11 10
3 5 4 5 5 5 5 4 5 4 4 4 5
4 8 8 7 8 8 8 8 8 7 7 8 8
5 9 9 9 9 9 9 8 8 8 8 8 9
6 3 3 3 3 3 3 3 3 2 3 3 3
str(timss.can)
'data.frame': 820 obs. of 84 variables:
$ IDCNTRY : Factor w/ 1 level " 90": 1 1 1 1 1 1 1 1 1 1 ...
$ IDSCH : Factor w/ 448 levels " 101"," 102",..: 1 1 2 3 3 4 5 5 6 7 ...
$ IDSCHOOL: Factor w/ 448 levels " 189"," 190",..: 307 307 308 309 309 310 311 311 312 313 ...
$ IDSTUD : Factor w/ 820 levels " 1890110"," 1900413",..: 520 521 522 523 524 525 526 527 528 529 ...
$ BC4GSBED: num 0 0 2 2 2 2 2 2 0 3 ...
$ BC4GSBEA: num 3 3 0 1 1 0 1 1 3 0 ...
$ BC4GCHTS: num 2 2 2 1 1 2 1 1 1 2 ...
$ BC4GCHES: num 1 1 1 1 1 2 1 1 1 1 ...
$ BC4GCHPS: num 1 1 1 2 2 2 2 2 2 2 ...
$ BCDGAS : num 1 1 0 1 1 2 2 2 2 2 ...
$ BCDSRMI : num 0 0 0 0 0 0 0 0 1 1 ...
$ BCDGPPSC: num 1 1 1 1 1 2 1 1 1 1 ...
$ ITLANG : Factor w/ 2 levels " 0"," 1": 2 2 2 2 2 2 2 2 2 2 ...
$ ITSEX : Factor w/ 2 levels " 0"," 1": 2 2 2 1 2 2 2 2 2 2 ...
$ BS4GBOOK: num 3 0 4 0 1 1 2 3 2 2 ...
$ BS4MHCAL: num 0 0 0 0 0 2 1 1 1 2 ...
$ BSDGEDUP: num 2 2 0 1 2 3 1 0 0 0 ...
$ BSDGCAVL: num 0 4 0 2 0 0 0 0 0 1 ...
$ BSDMPATM: num 2 2 1 0 0 1 2 1 1 1 ...
$ BSDMSVM : num 0 0 0 0 0 0 0 0 1 0 ...
$ BSDMSCM : num 2 2 0 1 1 1 1 1 0 0 ...
$ BSDGPBSS: num 2 1 0 0 1 1 2 0 0 1 ...
$ BC4GTENR: num 146 146 210 262 262 245 305 305 NA 430 ...
$ BCDGHW : num 28 28 25 26 26 26 26 26 NA 26 ...
$ BSDAGE : num 14 14 13 14 14 13 14 14 14 14 ...
$ BSMMAT05: num 563 500 548 475 521 482 474 545 540 460 ...
$ m1 : num 1 1 1 1 1 1 NA 1 1 1 ...
$ m2 : num 1 1 1 1 1 0 NA 1 1 1 ...
$ m3 : num 1 1 0 1 0 1 NA 0 1 0 ...
$ m4 : num 0 0 0 0 1 0 NA 0 1 0 ...
$ m5 : num 0 0 0 0 0 0 NA 0 0 1 ...
$ m6 : num 1 1 1 1 0 1 NA 0 0 1 ...
$ m7 : num 0 1 1 0 1 0 NA 1 1 0 ...
$ m8 : num 0 0 0 0 0 0 NA 0 1 0 ...
$ m9 : num 0 0 0 0 1 0 NA 1 1 0 ...
$ m10 : num 0 0 2 0 0 0 NA 1 0 0 ...
$ m11 : num 2 0 2 0 1 0 NA 0 2 0 ...
$ m12 : num 2 0 2 0 2 0 NA 0 0 0 ...
$ m13 : num 0 0 1 1 0 1 NA 0 0 0 ...
$ m14 : num 1 1 1 1 1 1 NA 1 1 0 ...
$ m15 : num 1 0 1 0 1 0 NA 1 0 0 ...
$ m16 : num 1 1 1 0 0 1 NA 1 0 1 ...
$ m17 : num 1 1 0 1 1 0 NA 0 0 1 ...
$ m18 : num 1 0 1 0 1 0 NA 1 1 1 ...
$ m19 : num 1 0 0 1 1 0 NA 1 1 0 ...
$ m20 : num 1 0 1 0 1 1 NA 1 1 0 ...
$ m21 : num 1 0 1 0 1 1 NA 1 1 0 ...
$ m22 : num 0 0 1 0 0 1 NA 0 0 0 ...
$ m23 : num 1 0 1 0 1 0 NA 0 1 1 ...
$ m24 : num 1 0 1 0 0 0 NA 0 0 0 ...
$ m25 : num 1 0 0 0 0 0 NA 1 1 0 ...
$ m26 : num 1 1 1 1 1 1 NA 1 1 0 ...
$ m27 : num 1 1 1 1 1 1 NA 1 1 1 ...
$ m28 : num 1 1 1 1 1 1 NA 0 1 0 ...
$ m29 : num 2 0 2 2 1 0 NA 0 0 0 ...
$ tot1 : num 23 10 24 12 19 11 NA 14 18 8 ...
$ tot2 : num 23 10 24 12 19 12 NA 14 18 8 ...
$ tot3 : num 23 10 25 12 20 11 NA 15 18 9 ...
$ tot4 : num 24 11 25 13 19 12 NA 15 18 9 ...
$ tot5 : num 24 11 25 13 20 12 NA 15 19 8 ...
$ tot6 : num 23 10 24 12 20 11 NA 15 19 8 ...
$ tot7 : num 24 10 24 13 19 12 NA 14 18 9 ...
$ tot8 : num 24 11 25 13 20 12 NA 15 18 9 ...
$ tot9 : num 24 11 25 13 19 12 NA 14 18 9 ...
$ tot10 : num 24 11 23 13 20 12 NA 14 19 9 ...
$ tot11 : num 22 11 23 13 19 12 NA 15 17 9 ...
$ tot12 : num 22 11 23 13 18 12 NA 15 19 9 ...
$ tot13 : num 24 11 24 12 20 11 NA 15 19 9 ...
$ tot14 : num 23 10 24 12 19 11 NA 14 18 9 ...
$ tot15 : num 23 11 24 13 19 12 NA 14 19 9 ...
$ tot16 : num 23 10 24 13 20 11 NA 14 19 8 ...
$ tot17 : num 23 10 25 12 19 12 NA 15 19 8 ...
$ tot18 : num 23 11 24 13 19 12 NA 14 18 8 ...
$ tot19 : num 23 11 25 12 19 12 NA 14 18 9 ...
$ tot20 : num 23 11 24 13 19 11 NA 14 18 9 ...
$ tot21 : num 23 11 24 13 19 11 NA 14 18 9 ...
$ tot22 : num 24 11 24 13 20 11 NA 15 19 9 ...
$ tot23 : num 23 11 24 13 19 12 NA 15 18 8 ...
$ tot24 : num 23 11 24 13 20 12 NA 15 19 9 ...
$ tot25 : num 23 11 25 13 20 12 NA 14 18 9 ...
$ tot26 : num 23 10 24 12 19 11 NA 14 18 9 ...
$ tot27 : num 23 10 24 12 19 11 NA 14 18 8 ...
$ tot28 : num 23 10 24 12 19 11 NA 15 18 9 ...
$ tot29 : num 22 11 23 11 19 12 NA 15 19 9 ...
- attr(*, "codepage")= int 65001
head(timss.can)
IDCNTRY IDSCH IDSCHOOL IDSTUD BC4GSBED BC4GSBEA BC4GCHTS BC4GCHES
1 90 101 3141 31410110 0 3 2 1
2 90 101 3141 31410124 0 3 2 1
3 90 102 3142 31420107 2 0 2 1
4 90 103 3143 31430213 2 1 1 1
5 90 103 3143 31430227 2 1 1 1
6 90 104 3144 31440217 2 0 2 2
BC4GCHPS BCDGAS BCDSRMI BCDGPPSC ITLANG ITSEX BS4GBOOK BS4MHCAL
1 1 1 0 1 1 1 3 0
2 1 1 0 1 1 1 0 0
3 1 0 0 1 1 1 4 0
4 2 1 0 1 1 0 0 0
5 2 1 0 1 1 1 1 0
6 2 2 0 2 1 1 1 2
BSDGEDUP BSDGCAVL BSDMPATM BSDMSVM BSDMSCM BSDGPBSS BC4GTENR BCDGHW
1 2 0 2 0 2 2 146 28
2 2 4 2 0 2 1 146 28
3 0 0 1 0 0 0 210 25
4 1 2 0 0 1 0 262 26
5 2 0 0 0 1 1 262 26
6 3 0 1 0 1 1 245 26
BSDAGE BSMMAT05 m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14 m15 m16
1 14 563 1 1 1 0 0 1 0 0 0 0 2 2 0 1 1 1
2 14 500 1 1 1 0 0 1 1 0 0 0 0 0 0 1 0 1
3 13 548 1 1 0 0 0 1 1 0 0 2 2 2 1 1 1 1
4 14 475 1 1 1 0 0 1 0 0 0 0 0 0 1 1 0 0
5 14 521 1 1 0 1 0 0 1 0 1 0 1 2 0 1 1 0
6 13 482 1 0 1 0 0 1 0 0 0 0 0 0 1 1 0 1
m17 m18 m19 m20 m21 m22 m23 m24 m25 m26 m27 m28 m29 tot1 tot2 tot3 tot4
1 1 1 1 1 1 0 1 1 1 1 1 1 2 23 23 23 24
2 1 0 0 0 0 0 0 0 0 1 1 1 0 10 10 10 11
3 0 1 0 1 1 1 1 1 0 1 1 1 2 24 24 25 25
4 1 0 1 0 0 0 0 0 0 1 1 1 2 12 12 12 13
5 1 1 1 1 1 0 1 0 0 1 1 1 1 19 19 20 19
6 0 0 0 1 1 1 0 0 0 1 1 1 0 11 12 11 12
tot5 tot6 tot7 tot8 tot9 tot10 tot11 tot12 tot13 tot14 tot15 tot16 tot17
1 24 23 24 24 24 24 22 22 24 23 23 23 23
2 11 10 10 11 11 11 11 11 11 10 11 10 10
3 25 24 24 25 25 23 23 23 24 24 24 24 25
4 13 12 13 13 13 13 13 13 12 12 13 13 12
5 20 20 19 20 19 20 19 18 20 19 19 20 19
6 12 11 12 12 12 12 12 12 11 11 12 11 12
tot18 tot19 tot20 tot21 tot22 tot23 tot24 tot25 tot26 tot27 tot28 tot29
1 23 23 23 23 24 23 23 23 23 23 23 22
2 11 11 11 11 11 11 11 11 10 10 10 11
3 24 25 24 24 24 24 24 25 24 24 24 23
4 13 12 13 13 13 13 13 13 12 12 12 11
5 19 19 19 19 20 19 20 20 19 19 19 19
6 12 12 11 11 11 12 12 12 11 11 11 12
#########################################
#stratification propensity score matching
######################################
# stage one: Obtain prpns
# 1-level model
xpsmod1<-glm(ITLANG~ ITSEX+BS4GBOOK+BS4MHCAL+BSDGEDUP+BSDGCAVL+BSDMPATM+BSDMSVM+BSDMSCM+BSDGPBSS,
data=timss, family=binomial)
summary(xpsmod1)
Call:
glm(formula = ITLANG ~ ITSEX + BS4GBOOK + BS4MHCAL + BSDGEDUP +
BSDGCAVL + BSDMPATM + BSDMSVM + BSDMSCM + BSDGPBSS, family = binomial,
data = timss)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.917 -0.904 -0.528 1.060 2.557
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.0940 0.1431 -0.66 0.5113
ITSEX 1 0.1366 0.0956 1.43 0.1530
BS4GBOOK 0.2770 0.0396 7.00 2.5e-12 ***
BS4MHCAL -0.4796 0.0424 -11.32 < 2e-16 ***
BSDGEDUP -0.2050 0.0416 -4.92 8.6e-07 ***
BSDGCAVL -0.3780 0.0587 -6.44 1.2e-10 ***
BSDMPATM 0.3529 0.0642 5.50 3.8e-08 ***
BSDMSVM 0.1382 0.0994 1.39 0.1646
BSDMSCM -0.1317 0.0755 -1.74 0.0812 .
BSDGPBSS -0.1941 0.0667 -2.91 0.0036 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 3115.9 on 2353 degrees of freedom
Residual deviance: 2674.6 on 2344 degrees of freedom
AIC: 2695
Number of Fisher Scoring iterations: 4
ps.all.1<-fitted(xpsmod1,response=T)
densityplot(ps.all.1)
grp1<-cut(ps.all.1,breaks=c(0,quantile(ps.all.1, c(.1,.2,.3,.4,.5,.6,.7,.8,.9)),1))
levels(grp1)
[1] "(0,0.112]" "(0.112,0.172]" "(0.172,0.237]" "(0.237,0.301]"
[5] "(0.301,0.363]" "(0.363,0.429]" "(0.429,0.496]" "(0.496,0.566]"
[9] "(0.566,0.664]" "(0.664,1]"
d1<-subset(timss,grp1==levels(grp1)[1]);head(d1)
IDCNTRY IDSCH IDSCHOOL IDSTUD BC4GSBED BC4GSBEA BC4GCHTS BC4GCHES
2 48 1 1 10126 1 3 1 1
3 48 1 1 10314 1 3 1 1
15 48 4 4 40210 1 0 2 2
18 48 4 4 40330 1 0 2 2
23 48 6 6 60108 0 2 0 1
25 48 6 6 60220 0 2 0 1
BC4GCHPS BCDGAS BCDSRMI BCDGPPSC ITLANG ITSEX BS4GBOOK BS4MHCAL
2 2 2 1 1 0 0 1 3
3 2 2 1 1 0 0 1 3
15 3 2 1 2 0 0 0 3
18 3 2 1 2 0 0 1 3
23 1 1 1 0 0 1 1 2
25 1 1 1 0 0 1 0 2
BSDGEDUP BSDGCAVL BSDMPATM BSDMSVM BSDMSCM BSDGPBSS BC4GTENR BCDGHW
2 2 1 0 0 1 2 570 25
3 2 4 0 0 0 1 570 25
15 2 4 0 0 1 2 277 25
18 2 1 0 0 1 1 277 25
23 3 4 0 0 0 2 NA 38
25 3 2 0 0 1 0 NA 38
BSDAGE BSMMAT05 m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14 m15 m16
2 13 382 1 0 1 0 0 1 0 0 0 0 0 0 0 1 1 1
3 14 352 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15 14 399 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1
18 15 259 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1
23 14 433 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0
25 18 343 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0
m17 m18 m19 m20 m21 m22 m23 m24 m25 m26 m27 m28 m29 tot1 tot2 tot3 tot4
2 0 0 0 0 0 0 0 1 0 1 1 1 2 11 12 11 12
3 0 0 1 0 0 0 0 1 0 1 1 1 0 5 5 5 5
15 0 0 1 1 0 0 0 1 0 0 1 0 0 6 7 7 7
18 0 0 1 0 0 0 0 0 1 0 0 0 0 5 5 5 4
23 0 1 1 1 0 0 0 1 1 1 1 1 0 11 11 10 11
25 0 1 1 0 0 1 0 0 1 1 1 0 0 9 8 9 9
tot5 tot6 tot7 tot8 tot9 tot10 tot11 tot12 tot13 tot14 tot15 tot16
2 12 11 12 12 12 12 12 12 12 11 11 11
3 5 5 5 5 5 5 5 5 5 5 5 5
15 7 7 6 7 7 7 7 7 7 7 7 6
18 5 4 5 5 5 5 5 5 5 5 5 4
23 11 10 11 11 11 11 11 11 11 10 11 11
25 9 8 9 9 9 9 9 9 9 9 8 9
tot17 tot18 tot19 tot20 tot21 tot22 tot23 tot24 tot25 tot26 tot27 tot28
2 12 12 12 12 12 12 12 11 12 11 11 11
3 5 5 4 5 5 5 5 4 5 4 4 4
15 7 7 6 6 7 7 7 6 7 7 6 7
18 5 5 4 5 5 5 5 5 4 5 5 5
23 11 10 10 10 11 11 11 10 10 10 10 10
25 9 8 8 9 9 8 9 9 8 8 8 9
tot29
2 10
3 5
15 7
18 5
23 11
25 9
d2<-subset(timss,grp1==levels(grp1)[2])
d3<-subset(timss,grp1==levels(grp1)[3])
d4<-subset(timss,grp1==levels(grp1)[4])
d5<-subset(timss,grp1==levels(grp1)[5])
d6<-subset(timss,grp1==levels(grp1)[6])
d7<-subset(timss,grp1==levels(grp1)[7])
d8<-subset(timss,grp1==levels(grp1)[8])
d9<-subset(timss,grp1==levels(grp1)[9])
d10<-subset(timss,grp1==levels(grp1)[10])
#nrow(data1)
table(d1$ITLANG,d1$IDCNTRY)
48 90 422 634 818
0 62 1 14 47 104
1 0 2 0 5 1
table(d2$ITLANG,d2$IDCNTRY)
48 90 422 634 818
0 54 2 25 45 76
1 2 9 13 7 2
table(d3$ITLANG,d3$IDCNTRY)
48 90 422 634 818
0 57 11 16 62 43
1 2 23 5 15 1
table(d4$ITLANG,d4$IDCNTRY)
48 90 422 634 818
0 36 18 20 53 47
1 1 27 13 20 4
table(d5$ITLANG,d5$IDCNTRY)
48 90 422 634 818
0 30 19 23 39 29
1 4 47 10 21 10
table(d6$ITLANG,d6$IDCNTRY)
48 90 422 634 818
0 19 45 17 27 30
1 1 54 10 23 9
table(d7$ITLANG,d7$IDCNTRY)
48 90 422 634 818
0 10 37 25 28 26
1 4 66 10 19 11
table(d8$ITLANG,d8$IDCNTRY)
48 90 422 634 818
0 4 48 13 21 15
1 2 95 5 15 17
table(d9$ITLANG,d9$IDCNTRY)
48 90 422 634 818
0 3 55 9 21 10
1 6 95 4 17 16
table(d10$ITLANG,d10$IDCNTRY)
48 90 422 634 818
0 0 44 9 14 7
1 5 122 8 18 8
# 2-level model: can use lmer or glmer packages
xpsmod<-glmer(ITLANG~ BC4GSBED+BC4GSBEA+BC4GCHTS+BC4GCHES+BC4GCHPS+BCDGAS+BCDSRMI+BCDGPPSC
+ ITSEX+BS4GBOOK+BS4MHCAL+BSDGEDUP+BSDGCAVL+BSDMPATM+BSDMSVM+BSDMSCM+BSDGPBSS
+ (1|IDSCH), data=timss, family=binomial)
Warning: failure to converge in 10000 evaluations
# ONE COUNTRY DATA
xpsmod<-glmer(ITLANG~ BC4GSBED+BC4GSBEA+BC4GCHTS+BC4GCHES+BC4GCHPS+BCDGAS+BCDSRMI+BCDGPPSC
+ ITSEX+BS4GBOOK+BS4MHCAL+BSDGEDUP+BSDGCAVL+BSDMPATM+BSDMSVM+BSDMSCM+BSDGPBSS
+ (1|IDSCH), data=timss.can, family=binomial)
summary(xpsmod)
Generalized linear mixed model fit by maximum likelihood ['glmerMod']
Family: binomial ( logit )
Formula: ITLANG ~ BC4GSBED + BC4GSBEA + BC4GCHTS + BC4GCHES + BC4GCHPS + BCDGAS + BCDSRMI + BCDGPPSC + ITSEX + BS4GBOOK + BS4MHCAL + BSDGEDUP + BSDGCAVL + BSDMPATM + BSDMSVM + BSDMSCM + BSDGPBSS + (1 | IDSCH)
Data: timss.can
AIC BIC logLik deviance
631.8 721.2 -296.9 593.8
Random effects:
Groups Name Variance Std.Dev.
IDSCH (Intercept) 45.2 6.72
Number of obs: 820, groups: IDSCH, 448
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.9396 2.0565 2.89 0.00387 **
BC4GSBED -1.0197 0.5410 -1.88 0.05947 .
BC4GSBEA -0.8135 0.4358 -1.87 0.06198 .
BC4GCHTS -0.6793 0.8264 -0.82 0.41111
BC4GCHES 1.1376 0.8444 1.35 0.17794
BC4GCHPS -2.8734 0.8539 -3.36 0.00077 ***
BCDGAS 1.5311 0.8216 1.86 0.06237 .
BCDSRMI 0.8321 0.8842 0.94 0.34664
BCDGPPSC 0.3855 1.2562 0.31 0.75895
ITSEX 1 0.4882 0.6424 0.76 0.44721
BS4GBOOK 0.5625 0.3045 1.85 0.06466 .
BS4MHCAL 0.7606 0.4480 1.70 0.08959 .
BSDGEDUP 0.2799 0.3749 0.75 0.45532
BSDGCAVL -0.2804 0.5507 -0.51 0.61060
BSDMPATM 0.2917 0.4350 0.67 0.50254
BSDMSVM -0.2641 0.6872 -0.38 0.70071
BSDMSCM -0.4011 0.5206 -0.77 0.44105
BSDGPBSS 0.0904 0.4937 0.18 0.85474
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) BC4GSBED BC4GSBEA BC4GCHT BC4GCHE BC4GCHP BCDGAS BCDSRM
BC4GSBED -0.164
BC4GSBEA -0.508 0.252
BC4GCHTS -0.199 0.018 -0.052
BC4GCHES -0.151 0.046 0.070 -0.171
BC4GCHPS -0.400 -0.222 0.193 -0.032 -0.108
BCDGAS -0.167 -0.219 0.004 -0.100 -0.032 -0.054
BCDSRMI -0.164 -0.036 -0.014 -0.067 0.084 -0.027 0.016
BCDGPPSC 0.068 0.011 0.047 -0.268 -0.332 -0.416 -0.092 -0.086
ITSEX1 -0.186 0.001 -0.074 -0.047 0.084 -0.031 0.008 0.034
BS4GBOOK -0.403 -0.021 -0.029 -0.026 0.057 0.017 0.023 0.066
BS4MHCAL -0.159 0.068 -0.057 -0.020 0.004 -0.069 0.012 0.017
BSDGEDUP -0.172 -0.021 0.106 0.041 0.004 -0.041 -0.030 -0.007
BSDGCAVL -0.055 -0.026 0.031 -0.032 0.004 0.007 -0.053 0.048
BSDMPATM -0.119 0.002 0.005 0.052 0.035 -0.052 -0.050 0.007
BSDMSVM -0.037 0.055 0.036 -0.047 -0.019 -0.017 0.024 0.074
BSDMSCM 0.024 -0.035 -0.031 -0.035 -0.043 0.063 -0.017 0.000
BSDGPBSS -0.130 0.061 -0.023 0.042 0.012 -0.010 -0.004 -0.036
BCDGPP ITSEX1 BS4GBO BS4MHC BSDGED BSDGCA BSDMPA BSDMSV BSDMSC
BC4GSBED
BC4GSBEA
BC4GCHTS
BC4GCHES
BC4GCHPS
BCDGAS
BCDSRMI
BCDGPPSC
ITSEX1 0.026
BS4GBOOK 0.078 0.091
BS4MHCAL 0.041 0.128 0.058
BSDGEDUP 0.011 0.034 0.227 -0.052
BSDGCAVL -0.021 -0.064 0.014 -0.007 -0.059
BSDMPATM -0.011 0.109 0.005 0.093 0.059 0.070
BSDMSVM 0.007 -0.101 0.003 0.037 -0.039 0.015 -0.173
BSDMSCM -0.015 -0.156 0.014 -0.092 -0.224 -0.042 -0.472 -0.015
BSDGPBSS -0.039 0.065 -0.051 0.060 -0.097 0.028 0.004 -0.048 -0.033
match_on.examples <- list()
#match_on.examples$ps1 <- match_on(xpsmod)
ps<-fitted(xpsmod,response=T);ps[1:10]
[1] 0.9984 0.9709 0.9905 0.9311 0.9899 0.9983 0.9930 0.9916 0.9942 0.9918
ps2<-predict(xpsmod);ps2[1:10]
[1] 6.407 3.508 4.649 2.603 4.585 6.373 4.958 4.768 5.152 4.798
densityplot(ps)
match_on.examples$ps <- match_on(ITLANG~ps, data=timss.can)
#fm1 <- fullmatch(match_on.examples$ps1, data = timss.can)
fm2 <- fullmatch(match_on.examples$ps, data = timss.can)
all.equal(names(fm2), row.names(timss.can))
[1] TRUE
timss.can.new<-cbind(timss.can, matches = fm2[row.names(timss.can)])
timss<-timss.can.new
###########################
##check distribution balance
############################
grp<-cut(ps,breaks=c(0,quantile(ps, c(.1,.2,.3,.4,.5,.6,.7,.8,.9)),1))
levels(grp)
[1] "(0,0.03636]" "(0.03636,0.07042]" "(0.07042,0.1227]"
[4] "(0.1227,0.9709]" "(0.9709,0.9878]" "(0.9878,0.9931]"
[7] "(0.9931,0.9967]" "(0.9967,0.9984]" "(0.9984,0.9996]"
[10] "(0.9996,1]"
data1<-subset(timss,grp==levels(grp)[1]);head(data1)
IDCNTRY IDSCH IDSCHOOL IDSTUD BC4GSBED BC4GSBEA BC4GCHTS
88 90 156 3199 31990704 1 2 2
89 90 156 3199 31990807 1 2 2
351 90 322 780 7800416 1 3 1
358 90 328 789 7890105 2 0 2
359 90 328 789 7890123 2 0 2
362 90 334 796 7960225 1 2 1
BC4GCHES BC4GCHPS BCDGAS BCDSRMI BCDGPPSC ITLANG ITSEX BS4GBOOK
88 1 3 0 0 1 0 1 3
89 1 3 0 0 1 0 0 3
351 1 2 0 1 1 0 1 2
358 2 4 1 0 2 0 0 0
359 2 4 1 0 2 0 1 1
362 1 2 1 0 1 0 0 0
BS4MHCAL BSDGEDUP BSDGCAVL BSDMPATM BSDMSVM BSDMSCM BSDGPBSS BC4GTENR
88 2 0 0 0 0 0 0 675
89 3 0 0 2 0 2 0 675
351 1 0 0 0 1 0 1 NA
358 1 2 1 1 0 1 2 105
359 0 0 0 0 0 0 1 105
362 2 0 1 1 0 2 1 419
BCDGHW BSDAGE BSMMAT05 m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14
88 37 13 570 1 1 1 1 1 1 1 0 1 2 2 0 1 1
89 37 13 488 1 1 1 0 0 0 1 1 0 0 1 0 0 1
351 25 14 558 0 0 1 0 1 1 0 1 1 0 0 0 1 1
358 25 15 429 0 1 0 1 0 0 1 0 0 0 0 0 0 1
359 25 14 525 1 1 1 1 1 0 0 1 1 2 0 0 1 1
362 25 14 421 0 0 0 0 0 0 0 0 0 0 0 0 0 1
m15 m16 m17 m18 m19 m20 m21 m22 m23 m24 m25 m26 m27 m28 m29 tot1 tot2
88 1 1 1 1 1 0 0 0 1 1 0 0 1 1 2 24 24
89 1 1 1 0 0 1 0 0 0 1 0 1 1 1 2 16 16
351 1 1 1 1 0 0 0 0 0 1 0 1 1 1 2 17 17
358 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 5 4
359 1 0 0 0 1 0 0 0 0 1 1 1 0 1 0 16 16
362 1 1 0 1 0 0 0 0 0 0 1 1 1 1 0 8 8
tot3 tot4 tot5 tot6 tot7 tot8 tot9 tot10 tot11 tot12 tot13 tot14 tot15
88 24 24 24 24 24 25 24 23 23 25 24 24 24
89 16 17 17 17 16 16 17 17 16 17 17 16 16
351 16 17 16 16 17 16 16 17 17 17 16 16 16
358 5 4 5 5 4 5 5 5 5 5 5 4 5
359 16 16 16 17 17 16 16 15 17 17 16 16 16
362 8 8 8 8 8 8 8 8 8 8 8 7 7
tot16 tot17 tot18 tot19 tot20 tot21 tot22 tot23 tot24 tot25 tot26
88 24 24 24 24 25 25 25 24 24 25 25
89 16 16 17 17 16 17 17 17 16 17 16
351 16 16 16 17 17 17 17 17 16 17 16
358 4 5 5 5 5 5 5 5 5 5 5
359 17 17 17 16 17 17 17 17 16 16 16
362 7 8 7 8 8 8 8 8 8 7 7
tot27 tot28 tot29 matches
88 24 24 23 1.497
89 16 16 15 1.106
351 16 16 15 1.436
358 5 5 5 1.321
359 17 16 17 1.351
362 7 7 8 1.335
data2<-subset(timss,grp==levels(grp)[2])
data3<-subset(timss,grp==levels(grp)[3])
data4<-subset(timss,grp==levels(grp)[4])
data5<-subset(timss,grp==levels(grp)[5])
data6<-subset(timss,grp==levels(grp)[6])
data7<-subset(timss,grp==levels(grp)[7])
data8<-subset(timss,grp==levels(grp)[8])
data9<-subset(timss,grp==levels(grp)[9])
data10<-subset(timss,grp==levels(grp)[10])
#nrow(data1)
t1<-table(data1$ITLANG,data1$IDCNTRY)
t2<-table(data2$ITLANG,data2$IDCNTRY)
t3<-table(data3$ITLANG,data3$IDCNTRY)
t4<-table(data4$ITLANG,data4$IDCNTRY)
t5<-table(data5$ITLANG,data5$IDCNTRY)
t6<-table(data6$ITLANG,data6$IDCNTRY)
t7<-table(data7$ITLANG,data7$IDCNTRY)
t8<-table(data8$ITLANG,data8$IDCNTRY)
t9<-table(data9$ITLANG,data9$IDCNTRY)
t10<-table(data10$ITLANG,data10$IDCNTRY)
cbind(t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
90 90 90 90 90 90 90
0 82 82 82 34 0 0 0
1 0 0 0 48 82 82 82
90 90 90
0 0 0 0
1 82 82 82
######################################################
# stage 2:logistic regression using stratified data
#str(xpsmod)
est1<-se1<-p1<-rep(NA,10)
est2<-se2<-p2<-rep(NA,10)
est3<-se3<-p3<-rep(NA,10)
for (i in 1:10) {
ft<-glm(m4~English+tot4+English*tot4+(1|IDSCH), data=timss3, subset=(grp==levels(grp)[i]), family=binomial)
ft<-summary(ft)
est1[i]<-coef(ft)[2];est2[i]<-coef(ft)[3];est3[i]<-coef(ft)[4]
se1[i]<-coef(ft)[2,2];se2[i]<-coef(ft)[3,2];se3[i]<-coef(ft)[4,2]
p1[i]<-coef(ft)[2,4];p2[i]<-coef(ft)[3,4];p3[i]<-coef(ft)[4,4]
}
Error: object 'timss3' not found
# parameter estimates and s.e. for 5 strata and for each predictor
data.frame(est1,se1,p1)
est1 se1 p1
1 NA NA NA
2 NA NA NA
3 NA NA NA
4 NA NA NA
5 NA NA NA
6 NA NA NA
7 NA NA NA
8 NA NA NA
9 NA NA NA
10 NA NA NA
data.frame(est2,se2,p2)
est2 se2 p2
1 NA NA NA
2 NA NA NA
3 NA NA NA
4 NA NA NA
5 NA NA NA
6 NA NA NA
7 NA NA NA
8 NA NA NA
9 NA NA NA
10 NA NA NA
data.frame(est3,se3,p3)
est3 se3 p3
1 NA NA NA
2 NA NA NA
3 NA NA NA
4 NA NA NA
5 NA NA NA
6 NA NA NA
7 NA NA NA
8 NA NA NA
9 NA NA NA
10 NA NA NA
# remove estimates from stratum-1 as the ratio is 3/525, almost no English test takers
est1r<-est1[2:10];est2r<-est2[2:10];est3r<-est3[2:10]
se1r<-se1[2:10];se2r<-se2[2:10];se3r<-se3[2:10]
# compute combined estimates based on inverse of variance (i.e., precision) weighting
est.comb1<-sum(est1r/se1r^2)/sum(1/se1r^2);se.comb1<-sqrt(1/sum(1/se1r^2))
est.comb2<-sum(est2r/se2r^2)/sum(1/se2r^2);se.comb2<-sqrt(1/sum(1/se2r^2))
est.comb3<-sum(est3r/se3r^2)/sum(1/se3r^2);se.comb3<-sqrt(1/sum(1/se3r^2))
c(est.comb1,est.comb2,est.comb3)
[1] NA NA NA
c(se.comb1,se.comb2,se.comb3)
[1] NA NA NA