capture log close
log using ch5_ordinal, replace  
* crea una bitácora con formato de texto plano

//  MODELOS PARA VARIABLES DEPENDIENTES ORDINALES

//  Regression Models for Categorical Dependent Variables - 2nd Edition
//  Chapter 5 - Models for Ordinal Outcomes
//  Long and Freese - 27Jul2005
//  Comentarios en español: Javier Aparicio - 3nov2008


version 9		  
* para que stata use sintaxis de version 9 (en caso de tener una versión posterior)
set scheme s2manual   // modo monocromático para gráficas
set more off		  // para que el output del DO file no haga pausa al correr


* Ejemplo de variable ordinal: 
* Actitudes frente a las madres que trabajan
* Var dep: Esta de acuerdo o no con la frase: 
* "Las mamas que trabajan pueden mantener relaciones calidas con sus hijos" 

sysuse ordwarm2, clear

desc warm yr89 male white age ed prst
sum warm yr89 male white age ed prst

tab warm
* Var. Dep.: ¿Está de acuerdo o no con la frase? 
* "Las mamas que trabajan pueden mantener relaciones calidas con sus hijos" 


//  Un modelo de regresión logístico ordinal 

ologit warm male white age ed prst, nolog	// muestra completa
* la opcion nolog omite el output de las iteraciones

ologit warm male white age ed prst if yr89==1, nolog  // muestra restringida


* comparando ORDINAL LOGIT con ORDINAL PROBIT 

ologit warm yr89 male white age ed prst, nolog
estimates store ologit  // guarda los estimadores ologit

oprobit warm yr89 male white age ed prst, nolog
estimates store oprobit  // guarda los estimadores oprobit

estimates table ologit oprobit, b(%9.3f) t label varwidth(30)
* crea una tabla con los modelos previos con tres decimales, t-stats y etiquetas de variables de hasta 30 caracteres


//  Pruebas de hipótesis para coeficientes individuales 

* Wald test
ologit warm male yr89 white age ed prst, nolog
test male   // test similar al z-score de la tabla

display "z*z=" -9.343*-9.343
return list
di "chi2=" r(chi2) "; sqrt(chi2)= " sqrt(r(chi2))

* LR (likelihood ratio) test  
ologit warm yr89 male white age ed prst, nolog // modelo "completo"
estimates store fullmod

* estimate constrained model
ologit warm yr89 white age ed prst, nolog  // un modelo restringido o anidado
estimates store nestedmod

lrtest fullmod nestedmod  		
* Compara el likelihood de ambos modelos
* Ho: el modelo anidado es igual de verosímil que el modelo completo

//  Pruebas de hipótesis para varios coeficientes 

ologit warm yr89 male white age ed prst, nolog

* Wald test
test age white male   
// Ho: las 3 variables son "conjuntamente insignificantes"

* LR test
ologit warm yr89 male white age ed prst, nolog  // modelo completo
estimates store fmodel

ologit warm yr89 ed prst, nolog  // modelo restringido
estimates store nmodel

lrtest fmodel nmodel
* Compara el likelihood de ambos modelos
* Ho: el modelo anidado es igual de verosímil que el modelo completo


//  Medidas de bondad de ajuste usando FITSTAT

ologit warm yr89 male white age ed prst, nolog
fitstat


//  Modelos OLOGIT y el supuesto de las "regresiones paralelas"
* Los modelos ologit asumen que cada variable X tiene el mismo
* efecto (positivo o negativo) en cada una de las categorías ordinales 
* de la variable dependiente: X afecta igual la prob(outcome1) 
* que prob(outomce2), etc., una vez que ajustamos por cada intercepto 
* o cutpoint.

* generamos datos para ilustrar el supuesto de regresiones paralelas
clear
set obs 41
gen n = 41-_n
replace n = (n-20)/2.5
gen p1 = exp((-1*n)-3.5)/(1+exp((-1*n)-3.5))
gen p2 = exp((-1*n))/(1+exp((-1*n)))
gen p3 = exp((-1*n)+2)/(1+exp((-1*n)+2))
gen x = (n+8)*(9.375)
label var p1 "Pr(y<=1 | x)"
label var p2 "Pr(y<=2 | x)"
label var p3 "Pr(y<=3 | x)"
graph twoway (line p1 p2 p3 x) , ylabel(0(.2)1) xlabel(0(50)150) ///
    ytitle("Pr(y<=m)") xtitle("x") ///
    ysize(2.6558) xsize(4.0413)
 
graph export 05parallel.emf, replace

* Una prueba LR aproximada del supuesto de regresiones paralelas 
* (proporcional odds assumption)

* nota: omodel es un comando extra de Stata: findit omodel

sysuse ordwarm2, clear
omodel logit warm yr89 male white age ed prst

* Una prueba de Wald del supuesto
ologit warm yr89 male white age ed prst, nolog
brant, detail

/*
brant performs a Brant test of the parallel regression assumption (also
called the proportional odds assumption) after ologit.  The test compares
slope coefficients of the J-1 binary logits implied by the ordered regression
model.  Stata reports both the results of an omnibus test for the entire
model and tests of the assumption for each of the independent variables in
the model.
*/

//  Identificando residuales y outliers con predict después de ologit

sysuse ordwarm2, clear
drop warmlt2 warmlt3 warmlt4
gen warmlt2 = (warm<2) if warm <. // omite recodificar los missing values
gen warmlt3 = (warm<3) if warm <.
gen warmlt4 = (warm<4) if warm <.

* regresión logit cuando warm < 2
logit warmlt2 yr89 male white age ed prst, nolog
predict rstd_lt2, rs  // guarda los "residuales estandarizados" de Pearson

* regresión logit cuando warm < 3
logit warmlt3 yr89 male white age ed prst, nolog
predict rstd_lt3, rs

* regresión logit cuando < 4
logit warmlt4 yr89 male white age ed prst, nolog
predict rstd_lt4, rs

* graficando los residuales cuando warm < 3

sort prst       // ordena los datos de acuerdo a prestigio
gen index = _n  // genera un índice de 1 al número de obs. 
                // (para el eje horizontal de la gráfica)

graph twoway scatter rstd_lt3 index, yline(0) ylabel(-4(2)4) ///
    xtitle("Observation Number")  xlabel(0(500)2293) ///
    ysize(2.6558) xsize(4.0413) msymbol(Oh)
graph export 05residplot.emf, replace


 

//  Calculando efectos estandarizados con LISTCOEF

sysuse ordwarm2, clear
ologit warm yr89 male white age ed prst, nolog
listcoef, std help

//  Graficando probabilidades predichas con predict

ologit warm yr89 male white age ed prst, nolog
predict SDwarm Dwarm Awarm SAwarm  
* guarda las predicciones de los 4 resultados posibles de warm

label var SDwarm "Pr(SD)"  // prob(Strongly Disagree)
label var Dwarm "Pr(D)"    // prob(Strongly Disagree)
label var Awarm "Pr(A)"    // prob(Strongly Disagree)
label var SAwarm "Pr(SA)"  // prob(Strongly Agree)

dotplot SDwarm Dwarm Awarm SAwarm, ylabel(0(.25).75) ///
        ysize(2.0124) xsize(3.039)
graph export 05pred.emf, replace

 

// Calculando probabilidades predichas para casos particulares con PRVALUE

ologit warm yr89 male white age ed prst, nolog

* Un obrero de 64 años encuestado en 1977…
prvalue, x(yr89=0 male=1 prst=20 age=64 ed=16) rest(mean)

* Una joven de 1989 muy educada y con un empleo muy prestigioso
prvalue, x(yr89=1 male=0 prst=80 age=30 ed=24) rest(mean) 

* Un individuo promedio encuestado en 1977
prvalue, x(yr89=0) rest(mean) brief
* la opción brief presenta un output resumido de prvalue

* Un individuo promedio encuestado en 1989
prvalue, x(yr89=1) rest(mean) brief


//  Tablas de probabilidades predichas con PRTAB

prtab yr89 male, novarlbl  
* Tabula las diferentes combinaciones de yr89 y male 
* en una tabla 2x2 (porque ambas son dummies)

* prtab no incluye intervalos de confianza, para obtenerlos
* hay que usar prvalue:

prvalue, x(yr89=0 male=1)
prvalue, x(yr89=1 male=0)


//  Cambios discretos con Intervalos de confianza

* Hombres vs. mujeres

prvalue , x(male=0 yr89=1) rest(mean) save // guarda la estimación
prvalue , x(male=1 yr89=1) rest(mean) dif  // calcula la diferencia con los  
                                           // resultados guardados

* Efecto de aumentar la edad en una desviación estándar 
* centrada en torno a la edad promedio (44.93)

display 44.935 - (.5*16.779)
display 44.935 + (.5*16.779)

prvalue, x(male=0 yr89=1 age=36.5455) rest(mean) save
prvalue, x(male=0 yr89=1 age=53.3245) rest(mean) dif

//  Graficando una serie de probabilidades predichas con PRGEN

sysuse ordwarm2, clear
ologit warm yr89 male white age ed prst, nolog

prgen age, from(20) to(80) generate(w89) x(male=0 yr89=1) ncases(13)

* Probabilidades predichas para edades de 20 a 80 años, usando 
* mujeres encuestadas en 1989 como línea de base.
* Las 13 predicciones se guardan en una serie de variables con el prefijo "w89"

desc w89*
label var w89p1 "SD"  // prob(outcome1)
label var w89p2 "D"   // prob(outcome2)
label var w89p3 "A"   // prob(outcome3)
label var w89p4 "SA"  // prob(outcome4)
label var w89s1 "SD"  
label var w89s2 "SD or D"  // probabilidad acumulada
label var w89s3 "SD, D or A"  // probabilidad acumulada

* Paso 1: Gráfica de probabilidades predichas

graph twoway connected w89p1 w89p2 w89p3 w89p4 w89x, ///
     title("Panel A: Predicted Probabilities") ///
     xtitle("Age") xlabel(20(10)80) ylabel(0(.25).50) ///
     yscale(noline) ylabel("") xline(44.93) ///
     ytitle("") name(tmp1, replace)

* Paso 2: Gráfica de probabilidades acumuladas

graph twoway connected w89s1 w89s2 w89s3 w89x, ///
    title("Panel B: Cumulative Probabilities") ///
    xtitle("Age") xlabel(20(10)80) ylabel(0(.25)1) xline(44.93) ///
    yscale(noline) ylabel("") name(tmp2, replace) ///
    ytitle("")

* Paso 3: combinar las gráficas guardadas en memoria
graph combine tmp1 tmp2, col(1) iscale(*.9) imargin(small) ///
    ysize(4.31) xsize(3.287)

* Paso 4: Guardar gráfica en formato "enhanced metafile"
graph export 05prgen.emf, replace


//  Calculando cambios en probabilidades predichas con PRCHANGE

ologit warm yr89 male white age ed prst, nolog

* Efectos marginales usando prchange (comando Spost)
prchange age, x(male=0 yr89=1) rest(mean) help

* Efectos marginales usando mfx (comando stata)
mfx compute, at(male=0 yr89=1) predict(outcome(1))

* Efectos marginales de cambios discretos con prchange
prchange male age prst, x(male=0 yr89=1) rest(mean) help

* Efecto marginal de un aumento de 10 años en edad
* (usando la opción delta de prchange)
quietly ologit warm yr89 male white age ed prst, nolog
prchange age, x(male=0 yr89=1) rest(mean) delta(10) help


//  Mostrar odds ratios con LISTCOEF

ologit warm yr89 male white age ed prst, nolog
listcoef male age, help
listcoef male age, percent help 
listcoef male, reverse help

log close