capture log close
log using irt, text replace
set more off
version 8.1

/*
Stata do-file for some of the analysis of Section 9.4 of
Skrondal, A. and Rabe-Hesketh, S (2004). Generalized
Latent Variable Modeling. Multilevel, Longitudinal
and Structural Equation Models. Chapman & Hall/CRC
*/

* read and prepare data
insheet using mislevy.dat, clear

* stack frequencies into wt2
gen i=_n
reshape long cw cb,i(i) j(male) string
replace i=_n
reshape long c, i(i) j(white) string
drop i

encode white, gen(w)
encode male, gen(m)
replace w=w-1
replace m=m-1

rename c wt2

* tot = number of people in group defined by w and m
egen tot = sum(wt2), by(w m)

* stack response into y and create variable item
gen patt=_n
reshape long y, i(patt) j(item)

* produce dummies for the item
qui tab item, gen(d)

** one-parameter logistic IRT model
gllamm y d1 d2 d3 d4, i(patt) l(logit) f(binom) weight(wt) nocons adapt

* item characeristic curve (x=standardized ability)
matrix list e(b)
twoway (function y=1/(1+exp(-[y]d1 -x*[patt1]_cons)), range(-2.5 2.5) clwidth(thick) )                  /*
*/     (function y=1/(1+exp(-[y]d2 -x*[patt1]_cons)), range(-2.5 2.5) clpatt(dot) clwidth(thick))       /*
*/     (function y=1/(1+exp(-[y]d3 -x*[patt1]_cons)), range(-2.5 2.5) clpatt(dash) clwidth(thick))      /*
*/     (function y=1/(1+exp(-[y]d4 -x*[patt1]_cons)), range(-2.5 2.5) clpatt(longdash) clwidth(thick)), /*
*/     legend( label(1 "Item 1") label(2 "Item 2") label(3 "Item 3") label(4 "Item 4") )  /*
*/     xtitle(Ability) ytitle(Probability of correct answer)

** two-parameter logisistic IRT model
eq load: d1-d4
gllamm y d1-d4, i(patt) eqs(load) l(logit) f(binom) weight(wt) nocons adapt nip(12) 

matrix list e(b)
twoway (function y=1/(1+exp(-[y]d1 -x*[pat1_1]d1)), range(-2.5 2.5) clwidth( thick))                              /*
*/     (function y=1/(1+exp(-[y]d2 -x*[pat1_1]d1*[pat1_1l]d2)), range(-2.5 2.5) clpatt(dot) clwidth(thick))      /*
*/     (function y=1/(1+exp(-[y]d3 -x*[pat1_1]d1*[pat1_1l]d3)), range(-2.5 2.5) clpatt(dash) clwidth(thick))      /*
*/     (function y=1/(1+exp(-[y]d4 -x*[pat1_1]d1*[pat1_1l]d4)), range(-2.5 2.5) clpatt(longdash) clwidth(thick)), /*
*/     legend( label(1 "Item 1") label(2 "Item 2") label(3 "Item 3") label(4 "Item 4"))            /*
*/     xtitle(Ability) ytitle(Probability of correct answer)

** two parameter logistic IRT with non-zero mean ability
gen cons=1
eq load: d1-d4
eq f1: cons
gllamm y d2-d4, i(patt) eqs(load) l(logit) f(binom) weight(wt) /*
  */ geqs(f1) nocons adapt nip(12)

* empirical Bayes predictions: ability scores
gllapred IRT, fac

** MIMIC model
gen f=1-m
gen b=1-w
gen b_f = b*f
eq f1: cons f b b_f
matrix a=e(b)
gllamm y d2-d4, i(patt) eqs(load) l(logit) f(binom) weight(wt) /*
  */ geqs(f1) from(a) nocons adapt nip(12)  

* empirical Bayes predictions: ability scores
gllapred MIMIC, fac

* look at ability scores for each response pattern
drop d1-d4
reshape wide y, i(patt) j(item)
sort y1-y4 b f
list y1-y4 b f IRTm1 MIMICm1, nolab clean


log close
exit