You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
415 lines
12 KiB
Plaintext
415 lines
12 KiB
Plaintext
8 months ago
|
program define irtpoly,eclass
|
||
|
|
version 11.0
|
||
|
|
syntax varlist(min=3 numeric) [if] [in] [,test Graph group(string) latent(string) REPlace Fixed(string) FIXEDVar(real -1) rsm rasch Last SASOUTput long Covariables(varlist) Covariablemean(varname) noCentered Project(string)]
|
||
|
|
preserve
|
||
|
|
|
||
|
|
|
||
|
|
capture mkdir "c:/data/irtpoly//`project'/"
|
||
|
|
if !_rc {
|
||
|
|
di in green "The directory c:/data/irtpoly//`project' has been created"
|
||
|
|
}
|
||
|
|
|
||
|
|
local dir="c:/data/irtpoly//`project'/"
|
||
|
|
local savegroup=1
|
||
|
|
if "`group'"=="" {
|
||
|
|
tempname group
|
||
|
|
local savegroup=0
|
||
|
|
}
|
||
|
|
local savelatent=1
|
||
|
|
if "`latent'"=="" {
|
||
|
|
tempname latent
|
||
|
|
local savelatent=0
|
||
|
|
}
|
||
|
|
tempvar order
|
||
|
|
gen `order'=_n
|
||
|
|
tempfile pcmsasfile
|
||
|
|
qui save `pcmsasfile'
|
||
|
|
qui count `if' `in'
|
||
|
|
local nbind=r(N)
|
||
|
|
tokenize `varlist'
|
||
|
|
local nbitems:word count `varlist'
|
||
|
|
local max=0
|
||
|
|
forvalues i=1/`nbitems' {
|
||
|
|
qui su ``i''
|
||
|
|
local max`i'=r(max)
|
||
|
|
if `max`i''>`max' {
|
||
|
|
local max=`max`i''
|
||
|
|
}
|
||
|
|
}
|
||
|
|
tempname freq
|
||
|
|
contract `varlist' `covariables' `covariablemean', freq(`freq')
|
||
|
|
qui sort `varlist'
|
||
|
|
qui outsheet `varlist' `covariables' `covariablemean' `freq' using "`dir'irtpoly_data.txt",replace
|
||
|
|
|
||
|
|
drop _all
|
||
|
|
if "`fixed'"!="" {
|
||
|
|
capture confirm matrix `fixed'
|
||
|
|
if _rc {
|
||
|
|
di as error "The {hi:`fixed'} matrix does not exist"
|
||
|
|
error 198
|
||
|
|
}
|
||
|
|
tempname matfix
|
||
|
|
matrix `matfix'=`fixed''
|
||
|
|
qui svmat `matfix'
|
||
|
|
qui rename `matfix'1 estimate
|
||
|
|
qui gen parameter=""
|
||
|
|
order parameter estimate
|
||
|
|
local l=1
|
||
|
|
forvalues j=1/`nbitems' {
|
||
|
|
forvalues m=1/`max`j'' {
|
||
|
|
if "`rasch'"=="" {
|
||
|
|
qui replace parameter="beta``j''_`m'" in `l'
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
qui replace parameter="beta``j''" in `l'
|
||
|
|
}
|
||
|
|
local ++l
|
||
|
|
}
|
||
|
|
}
|
||
|
|
qui outsheet using "`dir'/irtpoly_fixedparameters.txt",replace
|
||
|
|
|
||
|
|
}
|
||
|
|
drop _all
|
||
|
|
qui set obs 1000
|
||
|
|
qui generate str txt="%include 'C:\ado\macros SAS\anaqolv48.sas';" in 1
|
||
|
|
qui replace txt="%include 'C:\Documents and Settings\Jean-Benoit Hardouin\Mes documents\Boulot JB\Enseignement\ENSAI\2009-2010\macros\gammasymv1.sas';" in 2
|
||
|
|
qui replace txt="PROC IMPORT OUT=WORK.data DATAFILE='`dir'irtpoly_data.txt' DBMS=TAB REPLACE;GETNAMES=YES;DATAROW=2; RUN;" in 3
|
||
|
|
if "`fixed'"!="" {
|
||
|
|
qui replace txt="PROC IMPORT OUT=WORK.fixed DATAFILE='`dir'irtpoly_fixedparameters.txt' DBMS=TAB REPLACE;GETNAMES=YES;DATAROW=2; RUN;" in 4
|
||
|
|
}
|
||
|
|
local txt="%anaqol(DATASET=data,ITEMS=`varlist',DETAILS=yes,WEIGHT=`freq',MODEL="
|
||
|
|
if "`rsm'"==""&"`rasch'"=="" {
|
||
|
|
local txt `txt'pcm
|
||
|
|
}
|
||
|
|
else if "`rasch'"!="" {
|
||
|
|
local txt `txt'rasch, TEST=no
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
local txt `txt'rsm
|
||
|
|
}
|
||
|
|
if "`fixed'"!="" {
|
||
|
|
local txt `txt',FIXED=fixed
|
||
|
|
}
|
||
|
|
if `fixedvar'>0 {
|
||
|
|
local txt `txt',FIXEDVAR=`fixedvar'
|
||
|
|
}
|
||
|
|
if "`fixed'"!=""&"`fixedvar'"!="" {
|
||
|
|
local centered nocentered
|
||
|
|
}
|
||
|
|
if "`centered'"!="" {
|
||
|
|
local txt `txt',CENTERED=yes
|
||
|
|
}
|
||
|
|
if "`covariables'"!="" {
|
||
|
|
local txt `txt',COVARIABLES=`covariables'
|
||
|
|
}
|
||
|
|
if "`covariablemean'"!="" {
|
||
|
|
local txt `txt',COVARIABLEMEAN=`covariablemean'
|
||
|
|
}
|
||
|
|
local txt `txt');
|
||
|
|
qui replace txt="`txt'" in 10
|
||
|
|
qui replace txt="PROC EXPORT DATA= WORK.Out_parameters OUTFILE='`dir'irtpoly_parameters.txt' DBMS=TAB REPLACE;RUN;" in 11
|
||
|
|
qui replace txt="PROC EXPORT DATA= WORK.Out_latent OUTFILE='`dir'irtpoly_latent.txt' DBMS=TAB REPLACE;RUN;" in 12
|
||
|
|
qui replace txt="PROC EXPORT DATA= WORK.Out_rep OUTFILE='`dir'irtpoly_rep.txt' DBMS=TAB REPLACE;RUN;" in 13
|
||
|
|
qui replace txt="PROC EXPORT DATA= WORK.Out_fit OUTFILE='`dir'irtpoly_fit.txt' DBMS=TAB REPLACE;RUN;" in 14
|
||
|
|
|
||
|
|
qui outsheet txt using "`dir'irtpoly_pgmsas.txt", replace nonames noquote
|
||
|
|
|
||
|
|
if "`last'"=="" {
|
||
|
|
/*local date=c(current_date)
|
||
|
|
local jour=substr("`date'",1,2)
|
||
|
|
local mois=substr("`date'",4,3)
|
||
|
|
local an=substr("`date'",8,4)
|
||
|
|
if "`mois'"=="Jan" {local moisd 01}
|
||
|
|
if "`mois'"=="Feb" {local moisd 02}
|
||
|
|
if "`mois'"=="Mar" {local moisd 03}
|
||
|
|
if "`mois'"=="Apr" {local moisd 04}
|
||
|
|
if "`mois'"=="May" {local moisd 05}
|
||
|
|
if "`mois'"=="Jun" {local moisd 06}
|
||
|
|
if "`mois'"=="Jul" {local moisd 07}
|
||
|
|
if "`mois'"=="Aug" {local moisd 08}
|
||
|
|
if "`mois'"=="Sep" {local moisd 09}
|
||
|
|
if "`mois'"=="Oct" {local moisd 10}
|
||
|
|
if "`mois'"=="Nov" {local moisd 11}
|
||
|
|
if "`mois'"=="Dec" {local moisd 12}
|
||
|
|
di "`jour' `mois' `an' `moisd'"
|
||
|
|
shell "date" "01/01/2009"
|
||
|
|
*/
|
||
|
|
if "`long'"!=""{
|
||
|
|
local cmd winexec
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
local cmd shell
|
||
|
|
}
|
||
|
|
`cmd' "C:\Program Files\SAS\SAS 9.2\SASFoundation\9.2\sas.exe" "`dir'irtpoly_pgmsas.txt" -print "`dir'irtpoly_pgmsas.lst" -nolog
|
||
|
|
*shell "cmd.exe" "date `jour'/`moid'/`an'"
|
||
|
|
if "`long'"!="" {
|
||
|
|
exit
|
||
|
|
}
|
||
|
|
|
||
|
|
}
|
||
|
|
if "`sasoutput'"!="" {
|
||
|
|
view "`dir'irtpoly_pgmsas.lst"
|
||
|
|
}
|
||
|
|
*set trace on
|
||
|
|
|
||
|
|
*set trace on
|
||
|
|
drop _all
|
||
|
|
qui insheet using "`dir'irtpoly_fit.txt"
|
||
|
|
qui su value if descr=="-2 Log Likelihood"
|
||
|
|
local m2ll=r(mean)
|
||
|
|
local ll=-`m2ll'/2
|
||
|
|
qui su value if descr=="AIC (smaller is better)"
|
||
|
|
local aic=r(mean)
|
||
|
|
qui su value if descr=="BIC (smaller is better)"
|
||
|
|
local bic=r(mean)
|
||
|
|
|
||
|
|
drop _all
|
||
|
|
qui insheet using "`dir'irtpoly_parameters.txt"
|
||
|
|
tempname parameters separameters
|
||
|
|
qui su estimate if parameter=="var"
|
||
|
|
local variance=r(mean)
|
||
|
|
qui su standarderror if parameter=="var"
|
||
|
|
local sevariance=r(mean)
|
||
|
|
*set trace on
|
||
|
|
local nbcov:word count `covariables'
|
||
|
|
forvalues i=1/`nbcov' {
|
||
|
|
local cov`i':word `i' of `covariables'
|
||
|
|
qui su estimate if parameter=="beta`cov`i''"
|
||
|
|
local betacov`i'=r(mean)
|
||
|
|
qui su standarderror if parameter=="beta`cov`i''"
|
||
|
|
local secov`i'=r(mean)
|
||
|
|
}
|
||
|
|
*set trace off
|
||
|
|
*su
|
||
|
|
di in gr "Number of individuals: " in ye `nbind'
|
||
|
|
di in gr "Number of items: " in ye `nbitems'
|
||
|
|
di in gr "log-likelihood: " in ye %10.4f `ll'
|
||
|
|
di in gr "AIC: " in ye %10.4f `aic'
|
||
|
|
di in gr "BIC: " in ye %10.4f `bic'
|
||
|
|
di
|
||
|
|
di
|
||
|
|
if "`rsm'"=="" {
|
||
|
|
matrix `parameters'=J(`nbitems',`max',0)
|
||
|
|
matrix `separameters'=J(`nbitems',`max',0)
|
||
|
|
local l=1
|
||
|
|
|
||
|
|
forvalues i=1/`nbitems' {
|
||
|
|
forvalues j=1/`max`i'' {
|
||
|
|
if "`fixed'"=="" {
|
||
|
|
qui su estimate if parameter=="beta``i''_`j'"
|
||
|
|
matrix `parameters'[`i',`j']=r(mean)
|
||
|
|
qui su standarderror if parameter=="beta``i''_`j'"
|
||
|
|
matrix `separameters'[`i',`j']=r(mean)
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
matrix `parameters'[`i',`j']=`fixed'[1,`l']
|
||
|
|
}
|
||
|
|
local ++l
|
||
|
|
}
|
||
|
|
}
|
||
|
|
di in gr "{hline 52}"
|
||
|
|
di in gr "Items" _col(18) "Modality" _col(30) "Estimate" _col(39) "Standard error"
|
||
|
|
di in gr "{hline 52}"
|
||
|
|
forvalues j=1/`nbitems' {
|
||
|
|
di in gr "``j''" _c
|
||
|
|
forvalues m=1/`max`i'' {
|
||
|
|
di _col(25) in gr `m' _col(30) %8.4f in ye `parameters'[`j',`m'] _col(45) %8.4f in ye `separameters'[`j',`m']
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
matrix `parameters'=J(`=`nbitems'+`max'-1',1,0)
|
||
|
|
matrix `separameters'=J(`=`nbitems'+`max'-1',1,0)
|
||
|
|
local l=1
|
||
|
|
if "`fixed'"=="" {
|
||
|
|
forvalues i=1/`nbitems' {
|
||
|
|
qui su estimate if parameter=="beta``i''"
|
||
|
|
matrix `parameters'[`i',1]=r(mean)
|
||
|
|
qui su standarderror if parameter=="beta``i''"
|
||
|
|
matrix `separameters'[`i',1]=r(mean)
|
||
|
|
local ++l
|
||
|
|
}
|
||
|
|
forvalues l=`=`nbitems'+1'/`=`nbitems'+`max'-1' {
|
||
|
|
local m=`l'-`nbitems'+1
|
||
|
|
qui su estimate if parameter=="t`m'"
|
||
|
|
matrix `parameters'[`l',1]=r(mean)
|
||
|
|
local tau`m'=r(mean)
|
||
|
|
qui su standarderror if parameter=="t`m'"
|
||
|
|
matrix `separameters'[`l',1]=r(mean)
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
matrix `parameters'=`fixed'
|
||
|
|
}
|
||
|
|
di in gr "{hline 52}"
|
||
|
|
di in gr "Items" _col(30) "Estimate" _col(39) "Standard error"
|
||
|
|
di in gr "{hline 52}"
|
||
|
|
forvalues j=1/`nbitems' {
|
||
|
|
di in gr "``j''" _col(30) %8.4f in ye `parameters'[`j',1] _col(45) %8.4f in ye `separameters'[`j',1]
|
||
|
|
}
|
||
|
|
forvalues l=`=`nbitems'+1'/`=`nbitems'+`max'-1' {
|
||
|
|
di in gr "tau`=`l'-`nbitems''" _col(30) %8.4f in ye `parameters'[`l',1] _col(45) %8.4f in ye `separameters'[`l',1]
|
||
|
|
}
|
||
|
|
}
|
||
|
|
di in gr "{hline 52}"
|
||
|
|
di in gr "Variance" _col(30) %8.4f in ye `variance' _col(45) %8.4f in ye `sevariance'
|
||
|
|
di in gr "{hline 52}"
|
||
|
|
forvalues i=1/`nbcov' {
|
||
|
|
di in gr "`cov`i''" _col(30) %8.4f in ye `betacov`i'' _col(45) %8.4f in ye `secov`i''
|
||
|
|
}
|
||
|
|
if "`covariables'"!="" {
|
||
|
|
di in gr "{hline 52}"
|
||
|
|
}
|
||
|
|
*matrix list `parameters'
|
||
|
|
*fdsjklgvsjf
|
||
|
|
|
||
|
|
*set trace on
|
||
|
|
|
||
|
|
drop _all
|
||
|
|
qui insheet using "`dir'irtpoly_rep.txt"
|
||
|
|
qui sort anaqol_id
|
||
|
|
qui sort `varlist' `covariables' `covariablemean'
|
||
|
|
qui tempfile pcmsas
|
||
|
|
qui rename theta `latent'
|
||
|
|
qui rename stderrpred se`latent'
|
||
|
|
qui save `pcmsas',replace
|
||
|
|
qui use `pcmsasfile', clear
|
||
|
|
qui sort `varlist'
|
||
|
|
qui gen anaqol_id=_n
|
||
|
|
qui sort anaqol_id
|
||
|
|
qui sort `varlist' `covariables' `covariablemean'
|
||
|
|
|
||
|
|
/***********************************************
|
||
|
|
qui merge 1:1 anaqol_id using "`pcmsas'",nogen
|
||
|
|
***********************************************/
|
||
|
|
qui merge m:1 `varlist' `covariables' `covariablemean' using "`pcmsas'",nogen
|
||
|
|
|
||
|
|
*tempvar group
|
||
|
|
*set trace on
|
||
|
|
forvalues i=1/`nbcov' {
|
||
|
|
qui replace `latent'=`latent'+`betacov`i''*`cov`i''
|
||
|
|
}
|
||
|
|
*qui save `latent' using c:\latent.dta
|
||
|
|
qui gengroup `latent', det replace continuous newvariable(`group')
|
||
|
|
qui su `group'
|
||
|
|
local nbgroup=r(max)
|
||
|
|
forvalues g=1/`nbgroup' {
|
||
|
|
qui count if `group'==`g'
|
||
|
|
local group`g'=r(N)
|
||
|
|
}
|
||
|
|
forvalues i=1/`nbitems' {
|
||
|
|
*set trace on
|
||
|
|
tempname freq`i'
|
||
|
|
qui tab `group' ``i'',matcell(`freq`i'') row nofreq m
|
||
|
|
*matrix list `freq`i''
|
||
|
|
forvalues g=1/`nbgroup' {
|
||
|
|
qui count if `group'==`g'&``i''!=.
|
||
|
|
local freq`g'_`i'=r(N)
|
||
|
|
forvalues j=0/`max`i'' {
|
||
|
|
matrix `freq`i''[`g',`=`j'+1']=`freq`i''[`g',`=`j'+1']/`freq`g'_`i''
|
||
|
|
}
|
||
|
|
}
|
||
|
|
local D`i'=0
|
||
|
|
forvalues j=0/`max`i'' {
|
||
|
|
local D`i'_`j' exp(`j'*`latent'
|
||
|
|
forvalues l=1/`j' {
|
||
|
|
if "`rsm'"=="" {
|
||
|
|
local D`i'_`j' `D`i'_`j''-`parameters'[`i',`l']
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
local D`i'_`j' `D`i'_`j''-`parameters'[`i',1]
|
||
|
|
}
|
||
|
|
}
|
||
|
|
if "`rsm'"!="" {
|
||
|
|
forvalues m=2/`j' {
|
||
|
|
local D`i'_`j' `D`i'_`j''-`tau`m''
|
||
|
|
}
|
||
|
|
}
|
||
|
|
local D`i'_`j' `D`i'_`j'')
|
||
|
|
local D`i' `D`i''+`D`i'_`j''
|
||
|
|
}
|
||
|
|
}
|
||
|
|
tempvar theta2
|
||
|
|
qui gen `theta2'=0
|
||
|
|
forvalues g=1/`nbgroup' {
|
||
|
|
qui su `latent' if `group'==`g'
|
||
|
|
local thetag`g'=r(mean)
|
||
|
|
qui replace `theta2'=`thetag`g'' if `group'==`g'
|
||
|
|
}
|
||
|
|
local colors="blue red green gray pink purple"
|
||
|
|
|
||
|
|
*local chi2=0
|
||
|
|
forvalues i=1/`nbitems' {
|
||
|
|
local line`i'
|
||
|
|
local scatter`i'
|
||
|
|
tempvar propE``i'' propO``i''
|
||
|
|
qui gen `propE``i'''=0
|
||
|
|
qui gen `propO``i'''=0
|
||
|
|
forvalues j=0/`max`i'' {
|
||
|
|
local color:word `=`j'+1' of `colors'
|
||
|
|
tempvar propE``i''_`j' propO``i''_`j'
|
||
|
|
*matrix list `parameters'
|
||
|
|
*di "qui gen `propE``i''_`j''=`D`i'_`j''/(`D`i'')"
|
||
|
|
qui gen `propE``i''_`j''=`D`i'_`j''/(`D`i'')
|
||
|
|
*su `propE``i''_`j''
|
||
|
|
label variable `propE``i''_`j'' "Expected values / modality `j'"
|
||
|
|
local line`i' `line`i'' (line `propE``i''_`j'' `latent', lcolor(`color') lwidth(thick))
|
||
|
|
qui gen `propO``i''_`j''=0
|
||
|
|
forvalues g=1/`nbgroup' {
|
||
|
|
local tmp=`freq`i''[`g',`=`j'+1']
|
||
|
|
qui replace `propO``i''_`j''=`tmp' if `group'==`g'
|
||
|
|
}
|
||
|
|
label variable `propO``i''_`j'' "Observed values / modality `j'"
|
||
|
|
qui replace `propO``i'''=`propO``i'''+`j'*`propO``i''_`j''
|
||
|
|
qui replace `propE``i'''=`propE``i'''+`j'*`propE``i''_`j''
|
||
|
|
local scatter`i' `scatter`i'' (scatter `propO``i''_`j'' `theta2', mcolor(`color'))
|
||
|
|
}
|
||
|
|
qui sort `latent'
|
||
|
|
if "`graph'"!="" {
|
||
|
|
twoway `line`i'' `scatter`i'',name(``i'', replace)
|
||
|
|
}
|
||
|
|
label variable `propE``i''' "Expected values"
|
||
|
|
label variable `propO``i''' "Observed values"
|
||
|
|
if "`graph'"!="" {
|
||
|
|
twoway (line `propE``i''' `latent', lcolor(green) lwidth(thick)) (scatter `propO``i''' `theta2',mcolor(green)),name(``i''2, replace)
|
||
|
|
}
|
||
|
|
*set trace on
|
||
|
|
if "`test'"!="" {
|
||
|
|
local chi2=0
|
||
|
|
forvalues g=1/`nbgroup' {
|
||
|
|
qui ttest `propE``i'''=`propO``i''' if `group'==`g'
|
||
|
|
local t`g'=r(t)
|
||
|
|
qui count if `group'==`g'
|
||
|
|
local nb`g'=r(N)
|
||
|
|
di "local chi2=`chi2'+/*`nb`g''**/(`t`g'')^2"
|
||
|
|
local chi2=`chi2'+/*`nb`g''**/(`t`g'')^2
|
||
|
|
}
|
||
|
|
di "Chi-square statistics: " %8.4f `chi2'
|
||
|
|
local pchi2=chi2(`=`nbgroup'-1',`chi2')
|
||
|
|
di "p-values: " %8.4f `pchi2'
|
||
|
|
}
|
||
|
|
}
|
||
|
|
*set trace on
|
||
|
|
tempfile saveu
|
||
|
|
qui keep `order' `latent' se`latent' `group'
|
||
|
|
if `savegroup'==0 {
|
||
|
|
drop `group'
|
||
|
|
}
|
||
|
|
if `savelatent'==0 {
|
||
|
|
drop `latent'
|
||
|
|
drop se`latent'
|
||
|
|
}
|
||
|
|
sort `order'
|
||
|
|
qui save `saveu' ,replace
|
||
|
|
|
||
|
|
restore
|
||
|
|
qui gen `order'=_n
|
||
|
|
qui sort `order'
|
||
|
|
if "`replace'"!="" {
|
||
|
|
capture drop `group'
|
||
|
|
capture drop `latent'
|
||
|
|
capture drop se`latent'
|
||
|
|
}
|
||
|
|
qui merge 1:1 `order' using `saveu',nogen
|
||
|
|
end
|