*********** general data preparation *********
clear
set mem 100M
use charity, clear

*collapse data
gen one=1
collapse(sum) wt2=one, by(ta1-ta5)
gen id=_n

********** partial credit model ************
**data preparation
set more off
list in 1/2, clean
reshape long ta, i(id) j(item)
list in 1/10, clean
tab item, gen(d)
sort id item
drop if ta==.
sort id item
list in 1/8, clean
gen obs=_n
expand 4
sort id item obs
list in 1/8, clean

by obs, sort: gen x = _n-1
gen chosen = ta == x
tab item, gen(it)
list id-it2 in 1/8, clean

forvalues i=1/5{
 forvalues g=1/3{
  gen d`i'_`g' = -1*it`i'*(x>=`g')
 }
}
list d1_1-d2_3 in 1/8, clean

*1PL partial credit model
eq slope: x
gllamm x d1_1-d5_3, i(id) eqs(slope) link(mlogit) expand(obs chosen o) weight(wt) adapt trace nocons

*CPC curves for Item1
quietly egen N=max(id)
generate trait1 = (-4) + (id-1)*(4-(-4))/(N-1) /*generate a latent scale from -4 to 4*/
gllapred prob1, mu us(trait)

estimates store pcm

set scheme sj
twoway (line prob1 trait1 if x==0, sort)                          ///
       (line prob1 trait1 if x==1, sort lpatt("."))                 ///
       (line prob1 trait1 if x==2, sort lpatt("-"))               ///
       (line prob1 trait1 if x==3, sort lpatt("_")) if item==4,   ///
legend(order(1 "strongly agree" 2 "somewhat agree"                ///
               3 "somewhat disagree" 4 "strongly disagree"))

graph save CPC1

*2PL partial credit model
*forvalues i=1/5{
* gen x_it`i'=x*it`i'
*}
*sort id item x
*list x_it1 x_it2 in 1/8, clean

*eq load: x_it1-x_it5
*gllamm x d1_1-d5_3, i(id) eqs(load) link(mlogit) expand(obs chosen o) weight(wt) adapt trace nocons

************ rating scale model *************
*data preparation
list in 1/8, clean

gen step1 = -1*(x>=1)
gen step2 = -1*(x>=2)
gen step3 = -1*(x>=3)

foreach var of varlist it*{
 gen n`var' = -1*`var'*x
}

sort id item x
list nit1 nit2 step2 step3 in 1/8, clean

*1PL rating scale model
eq slope: x
gllamm x nit1-nit5 step2 step3, i(id) eqs(slope) link(mlogit) expand(obs chosen o) weight(wt) adapt trace nocons

estimates store rsm

*2PL rating scale model
*eq load: x_it1-x_it5
*gllamm x nit1-nit5 step2 step3, i(id) eqs(load) link(mlogit) expand(obs chosen o) weight(wt) adapt trace nocons

lrtest rsm pcm

*CPC curves for Item4
*quietly egen N=max(id)
replace trait1 = (-4) + (id-1)*(4-(-4))/(N-1)
gllapred prob2, mu us(trait)

twoway (line prob2 trait1 if x==0, sort)                       ///
       (line prob2 trait1 if x==1, sort lpatt("."))                ///
       (line prob2 trait1 if x==2, sort lpatt("-"))              ///
       (line prob2 trait1 if x==3, sort lpatt("_")) if item==4,  ///
legend(order(1 "strongly agree" 2 "somewhat agree"                ///
               3 "somewhat disagree" 4 "strongly disagree"))


graph save CPC2