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

//  MODELOS PARA VARIABLES DEPENDIENTES CATEGÓRICAS

//  Regression Models for Categorical Dependent Variables - 2nd Edition
//  Chapter 6 - Models for Nominal Outcomes for Case-based Data
//  Long and Freese - 11Aug2005
//  Comentarios en español: Javier Aparicio - 10nov2008


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

//  EL MODELO LOGISTICO MULTINOMIAL

* Construimos una variable dependiente con 2 o 3 categorías

sysuse nomocc2, clear  // Datos de la 1982 General Social Survey

tab occ
gen occ3 = occ
recode occ3 1/3=1 4=2 5=3  // recodificamos a sólo 3 categorías

label def locc3 1 Manual 2 WhiteCol 3 Prof  // define value labels
label val occ3 locc3                        // asigna value labels

tab occ occ3,miss
gen wc_man= occ3==2 if occ3~=3     // whitecollar vs. manual dummy
gen prof_man = occ3==3 if occ3~=2  // professional vs. manual dummy
gen prof_wc = occ3==3 if occ3~=1   // professional vs. whitecollar dummy

label def lwm 1 WhiteCol 0 Manual
label def lpm 1 Prof 0 Manual
label def lpw 1 Prof 0 WhiteCol
label val wc_man lwm
label val prof_man lpm
label val prof_wc lpw

tab occ3 wc_man, miss  // la opción missing tabula valores omitidos
tab occ3 prof_man, miss
tab occ3 prof_wc, miss

tab prof_man, miss

* Modelo logit para resultados binarios
tab prof_man, miss
logit prof_man ed, nolog

* estimando el modelo multinomial logit correspondiente
tab occ3, miss
mlogit occ3 ed, nolog // compara el coef de ed en logit vs. mlogit

* cambiar los valores de las categorías nominales no afecta la estimación
gen occ999 = occ3
recode occ999 1=-3 2=0 3=999  // recodificación
mlogit occ3 ed, nolog base(3)   // la 3a categoría es la línea de base
mlogit occ999 ed, nolog base(999) // los coef de ambos modelos son idénticos


//  Un modelo para tipos de ocupación

desc white ed exper
sum white ed exper
tab occ
mlogit occ white ed exper, base(5) nolog  // modelo para las 5 categorías

//  los efectos del modelo anterior pueden reexpresarse
//  para diferentes categorías base (pairwise transitions) con LISTCOEF

listcoef white, help  // efectos de la variable white
listcoef, pvalue(.05)  // lista sólo los coeficientes significativos al 5%
listcoef ed, pvalue(.05) gt nolabel
* lista sólo coeficientes positivos, significativos al 5% y sin value labels

//  Al igual que en logit, las variables exactamente correlacionadas con cierto
// resultado (perfect predictors) no pueden ser estimados por mlogit. Ejemplo:

sysuse ordwarm2, clear
gen dummyprest = (prst<20 & warm==1)  // una dummy perfectamente correlacionada con warm=1
tab dummyprest warm, miss

mlogit warm yr89 male white age ed dummyprest, nolog
* noten cómo la dummy no resulta significativa

drop if dummyprest==1
mlogit warm yr89 male white age ed, nolog

//  Pruebas de hipótesis para variables individuales
//  Noten que para cada variable ahora se estiman K-1 coeficientes
//  (K=número de categorías)
//  Los Wald y Likelihood-ratio tests son similares pero se calculan diferente

sysuse nomocc2, clear
quietly mlogit occ white ed exper, base(5) nolog

* Prueba de Wald para una sola variable con test
test white
test ed
test exper
* Pruebas de Wald con MLOGTEST
mlogtest, wald

* Pruebas de likelihood-ratio (LR tests) para una sola variable con lrtest
mlogit occ white ed exper, base(5) nolog
estimates store fmodel  // guarda el modelo completo

* test of white
quietly mlogit occ ed exper, base(5) nolog
estimates store nmodel_white   // guarda el modelo restringido (sin white)
lrtest fmodel nmodel_white     // compara ambos modelos

* test of ed
quietly mlogit occ white exper, base(5) nolog
estimates store nmodel_ed    // guarda el modelo restringido (sin ed)
lrtest fmodel nmodel_ed      // compara ambos modelos

* test of exper
quietly mlogit occ white ed, base(5) nolog
estimates store nmodel_exper
lrtest fmodel nmodel_exper

* Las tres pruebas anteriores se pueden hacer fácilmemente con MLOGTEST
quietly mlogit occ white ed exper, base(5) nolog
mlogtest, lr


* Pruebas de Wald para múltiples variables con test
mlogit occ white ed exper, base(5) nolog
test ed exper

* Prueba de significancia conjunta usando MLOGTEST
mlogit occ white ed exper, base(5) nolog
mlogtest, wald set(ed exper)

* Pruebas de likelihood-ratio (LR tests) para múltiples variables con lrtest
mlogit occ white ed exper, base(5) nolog
estimates store fmodel
mlogit occ white, base(5) nolog
estimates store nmodel  // modelo restringido (sin ed exper)
lrtest fmodel nmodel

* Prueba de significancia conjunta usando MLOGTEST
mlogit occ white ed exper, base(5) nolog
mlogtest, lr set(ed exper)

* ¿Es posible reducir el número de categorías?
* MLOGIT estima un conjunto de coeficientes para cada categoría
* Si los coeficientes entre dos categorías o más son suficientemente
* similares es plausible colapsar o combinar estas categorías

* Prueba de Wald para combinar categorías
mlogit occ white ed exper, base(5) nolog
test [Menial]  // Ho: ningún coeficientes para categoría Menial es significativo
test [Menial=Craft] // Ho: los coeficientes para Menial son iguales a los de Craft

mlogtest, combine  // Si no rechazo Ho las categorías pueden combinarse

* LR test para combinar categorías
mlogit occ white ed exper, base(5) nolog
mlogtest, lrcomb


//  MLOGIT asume "independencia de alternativas irrelevantes"
//  (independence of irrelevant alternatives, IIA)
//  Pero este supuesto se puede poner a prueba

* Prueba de IIA de Hausman
mlogit occ white ed exper, base(5) nolog
mlogtest, hausman base

* Y si cambiamos el grupo base los resultados son idénticos
quietly mlogit occ white ed exper, base(1) nolog
mlogtest, hausman base
quietly mlogit occ white ed exper, base(2) nolog
mlogtest, hausman base

* Prueba de IIA de Small-Hsiao
quietly mlogit occ white ed exper, base(5) nolog
set seed 339487731
mlogtest, smhsiao
set seed 8675309
mlogtest, smhsiao


//  Calculando probabilidades predichas usando predict

* estimamos el modelo y calculamos predicciones
mlogit occ white ed exper, base(5) nolog
predict ProbM ProbB ProbC ProbW ProbP
* se guardan las predicciones de los 5 resultados posibles de occ
desc Prob*
sum Prob*

* mlogit calcula diferentes coeficientes para cada grupo, mientras que
* ologit asume datos ordinales y calcula un sólo conjunto
* de coeficientes.
* Podemos usar predict para comparar mlogit vs. ologit

sysuse ordwarm2,clear

ologit warm yr89 male white age ed prst, nolog
predict SDologit Dologit Aologit SAologit  // predicciones ologit
label var Dologit "ologit-D"

mlogit warm yr89 male white age ed prst, nolog
predict SDmlogit Dmlogit Amlogit SAmlogit  // predicciones mlogit
label var Dmlogit "mlogit-D"

dotplot Dologit Dmlogit, ylabel(0(.25).75)      ///
     ysize(2.0254) xsize(3.0381)   // comparamos predicciones de "Disagree"
graph export 06omlogdot.emf, replace

* Otra forma de comparar predicciones de ologit vs mlogit
corr SDologit SDmlogit // predicciones para "StronglyDisagree"
corr Dologit Dmlogit   // predicciones para "Disagree"
corr Aologit Amlogit   // predicciones para "Agree"
corr SAologit SAmlogit // predicciones para "StronglyAagree"

* graficando las dos predicciones para Disagree
graph twoway scatter Dologit Dmlogit,       ///
    xtitle("mlogit: Pr(y=Disagree)")        ///
    ytitle("ologit: Pr(y=Disagree)")        ///
    ylabel(0(.25).75) xlabel(0(.25).75)  sort
graph export 06omloggraph.emf, replace


// Calculando probabilidades predichas para casos particulares con PRVALUE

sysuse nomocc2, clear
mlogit occ white ed exper, base(5) nolog
prvalue, x(white=0) rest(mean) save
prvalue, x(white=1) rest(mean) dif

//  Tablas de probabilidades predichas con PRTAB

label def lwhite 0 NonWhite 1 White
label val white lwhite
prtab ed white, novarlbl outcome(1)

//  Graficando una serie de probabilidades predichas con PRGEN

* Plot de probabilidades de dos grupos (white y no-white)
quietly mlogit occ white ed exper, base(5)
prgen ed, x(white=1) from(6) to(20) generate(wht) ncases(15)
desc wht*
* Probabilidades predichas para personas blancas con educación de 6 a 20 años
* y con experiencia promedio. Las 15 predicciones se guardan en una serie
* de variables con el prefijo "wht"

prgen ed, x(white=0) from(6) to(20) generate(nwht) ncases(15)
desc nwht*
* Probabilidades predichas para personas NO blancas con educación de 6 a 20
* años y con experiencia promedio. Las 15 predicciones se guardan en una serie
* de variables con el prefijo "nwht"

label var whtp1 "Whites"
label var nwhtp1 "Nonwhites"

graph twoway connected whtp1 nwhtp1 nwhtx,  ///
    xtitle("Years of Education")            ///
    ytitle("Pr(Menial Job)")                ///
    ylabel(0(.25).50) xlabel(6 8 12 16 20)  ///
    ysize(2.7051) xsize(4.0421)
graph export 06prmenial.emf, replace

* Gráfica de probabilidades predichas para todos los resultados
* del grupo de personas blancas
label var whts1 "Menial"
label var whts2 "Blue Collar"
label var whts3 "Craft"
label var whts4 "White Collar"
graph twoway connected whts1 whts2 whts3 whts4 whtx,    ///
    xtitle("Whites: Years of Education")                ///
    ytitle("Summed Probability")                        ///
    xlabel(6(2)20) ylabel(0(.25)1) ysize(2.6195) xsize(4.0421)
graph export 06prsum.emf, replace

* Misma gráfica anterior con sombreado de áreas (rarea)
gen zero = 0
gen one  = 1
graph twoway (rarea zero whts1 whtx, bc(gs1))           ///
    (rarea whts1 whts2 whtx, bc(gs4))                   ///
    (rarea whts2 whts3 whtx, bc(gs8))                   ///
    (rarea whts3 whts4 whtx, bc(gs11))                  ///
    (rarea whts4 one whtx, bc(gs14)),                   ///
    ytitle("Summed Probability")                        ///
    legend( order( 1 2 3 4 5)                           ///
    label( 1 "Menial")                                  ///
    label( 2 "Blue Collar") label( 3 "Craft")           ///
    label(4 "White Collar") label(5 "Professional"))    ///
    xtitle("Whites: Years of Education")                ///
    xlabel(6 8 12 16 20) ylabel(0(.25)1)                ///
    ysize(2.6195) xsize(4.0421)                         ///
    plotregion(margin(zero))
graph export 06prsumrevised.emf, replace

//  Calculando cambios en probabilidades predichas con PRCHANGE

mlogit occ white ed exper, base(5) nolog
* prchange
prchange
* mfx compute: ojo, este comando puede tardar algo de tiempo
mfx compute, predict(outcome(1))

** Hasta aquí llegamos en el curso. 
** Long y Freese cubren otros temas más abajo en el DO file.

log close 

STOP!!


//  Graficando cambios discretos en prob. predichas con MLOGPLOT
//  Nota: estos ejemplos pueden reproducirse con la interfaz MLOGVIEW

quietly mlogit occ white ed exper, base(5) nolog
mlogplot, std(0uu) p(.1) dc ntics(9)
mlogplot, std(0ss) p(.1) min(-.2) max(.4) dc ntics(7)


//  odds ratios using listcoef and mlogview (pg 224 RevEd)

* plot hypothetical values: relative to A
matrix mnlbeta = (-.693, .693, .347 \ .347, -.347, .693 )
matrix mnlsd = (1, 2, 4)
global mnlname = "x1 x2 x3"
global mnlcatnm = "B C A"
global mnldepnm "depvar"
mlogplot, matrix std(uuu) vars(x1 x2 x3) packed
graph export 06orploteg.emf, replace

* plot hypothetical values: relative to B
matrix mnlbeta = (-.693, .693, .347 \ .347, -.347, .693 )
matrix mnlsd = (1, 2, 4)
global mnlname = "x1 x2 x3"
global mnlcatnm = "B C A"
global mnldepnm "depvar"
* note: use base(1) since B is the first letter in mnlcatnm
mlogplot, matrix std(uuu) vars(x1 x2 x3) packed base(1)
graph export 06orplotegb.emf, replace

* plot odds ratios from occupational outcomes example
mlogit occ white ed exper,base(5) nolog
listcoef white, help
* or plot without significance levels
mlogplot white ed exper, std(0ss) b(5) p(1) min(-2.75) max(.55) or  ///
    packed ntics(7)
graph export 06orplotocc1.emf,replace

* or plot with significance levels
mlogplot white ed exper, prob(.1) std(0ss) b(5) min(-2.75) max(.55) ///
    or ntics(7)
graph export 06orplotocc2.emf, replace
* or plot with significance levels and discrete change
prchange
mlogplot white ed exper, prob(.1) std(0ss) b(5) min(-2.75) max(.55) ///
    dc or ntics(7)
graph export 06orplotocc3.emf, replace

//  using mlogplot (pg 230 RevEd)

* or plot using mlogplot
mlogit occ white ed exper, base(5) nolog
prchange
mlogplot  white ed exper, std(0ss) p(.1) min(-2.75) max(.55) or ntics(7)

* other examples using mlogplot (not in book)
mlogplot white ed exper, std(0ss) p(.1) min(-2.75) max(.55) or dc ntics(7)
mlogplot white ed exper, dc std(0ss) min(-.5) max(.5)
mlogplot white ed exper, or std(0ss) min(-2.5) max(.5)
mlogplot white ed exper, or dc std(0ss) min(-.5) max(.5)
mlogplot white ed exper, dc or std(0ss) min(-2.5) max(.5) p(.1)

//   plotting estimates from matrices with mlogplot (pg 231 RevEd)

* plot example data
matrix mnlbeta = (-.693, .693, .347 \ .347, -.347, .693 )
matrix mnlsd   = (1, 2, 4)
global mnlname = "x1 x2 x3"
global mnlcatnm = "B C A"
global mnldepnm "depvar"
mlogplot, matrix std(uuu) vars(x1 x2 x3) packed
graph export 06matrix1.emf, replace

//  combining coefficients from two groups

* estimate two sets of coefficients
sysuse nomocc2, clear
mlogit occ ed exper if white==1, base(5) nolog
mlogit occ ed exper if white==0, base(5) nolog
* save coefficients for ed to mnlbeta
matrix mnlbeta =                    ///
    (-.8307514, -.9225522, -.6876114, -.4196403 \   ///
    -.7012628, -.560695 , -.882502 , -.5311514 )
* transpose to make columns correspond to variables
matrix mnlbeta = mnlbeta'
* compute sd of ed and add to matrix
sum ed
matrix mnlsd = (2.946427,2.946427)
global mnlname = "White NonWhite"
global mnlcatnm = "Menial BlueCol Craft WhiteCol Prof"
mlogplot, vars(White NonWhite) packed               ///
    or matrix std(ss)                   ///
    note("Racial Differences in Effects of Education")
graph export 06matrix2.emf, replace

//  multinomial probit (pg  2Ed)

* compare mprobit coefficients to probit coefficients
sysuse binlfp2, clear
probit lfp k5 k618 age wc hc lwg inc, nolog
mprobit lfp k5 k618 age wc hc lwg inc, nolog base(0)
mprobit lfp k5 k618 age wc hc lwg inc, nolog base(0) probitparam

* compare mprobit and probit predictions
sysuse binlfp2, clear
probit lfp k5 k618 age wc hc lwg inc, nolog
predict p_probit1
mprobit lfp k5 k618 age wc hc lwg inc, nolog base(0)
predict p_mprobit0 p_mprobit1
mprobit lfp k5 k618 age wc hc lwg inc, nolog base(0) probitparam
predict p_mprobit0p p_mprobit1p
pwcorr p_probit1 p_mprobit1 p_mprobit0 p_mprobit1p p_mprobit0p

* compare mprobit and mlogit
sysuse nomocc2, clear
mprobit occ white ed exper, base(1)
predict mpp1 mpp2 mpp3 mpp4 mpp5
* with probit normalization
mprobit occ white ed exper, base(1) probitparam
mlogit occ white ed exper, base(1)
predict mlp1 mlp2 mlp3 mlp4 mlp5
corr *p1
gen prdif1 = abs(mpp1-mlp1)
gen prdif2 = abs(mpp2-mlp2)
gen prdif3 = abs(mpp3-mlp3)
gen prdif4 = abs(mpp4-mlp4)
gen prdif5 = abs(mpp5-mlp5)
sum prdif*

//  stereotype logit model (pg  2Ed)

sysuse ordwarm2, replace
tab warm
tab warm, nolabel

* mlogit for attitudes toward working mothers
mlogit warm yr89 male white age ed prst, nolog base(4)

* slm for attitudes toward working mothers
slogit warm yr89 male white age ed prst, nolog

* predictions for the estimation sample
predict slpr1 slpr2 slpr3 slpr4
list warm slpr1-slpr4 in 1/4
mlogit warm yr89 male white age ed prst, nolog base(4)
predict mlpr1 mlpr2 mlpr3 mlpr4
list warm slpr1 mlpr1 slpr2 mlpr2 slpr3 mlpr3 in 1/5
pwcorr slpr1 mlpr1
pwcorr slpr2 mlpr2
pwcorr slpr3 mlpr3
pwcorr slpr4 mlpr4

* using pr* commands
prvalue , x(male=0) rest(mean) save
prvalue , x(male=1) rest(mean) dif
prchange male, rest(mean)

* odds ratios
listcoef, expand

* distinguishability and the phi parameters
slogit warm yr89 male white age ed prst, nolog
est store slbase
* wald test that phi1=phi2
nlcom ([phi1_1]_b[_cons] - [phi1_2]_b[_cons])
* phi1=phi2
constraint define 1 [phi1_1]_cons=[phi1_2]_cons
slogit warm yr89 male white age ed prst, nolog constraint(1)
lrtest slbase .
* phi2=phi3
constraint define 2 [phi1_2]_cons=[phi1_3]_cons
slogit warm yr89 male white age ed prst, nolog constraint(2)
lrtest slbase .
* phi3=phi4
constraint define 2 [phi1_3]_cons=0
slogit warm yr89 male white age ed prst, nolog constraint(2)
lrtest slbase .
* phi1=phi2 and phi3=phi4
slogit warm yr89 male white age ed prst, nolog constraint(1 2)
lrtest slbase .

* ordinality and the 1-dimensional slogit model
gen warm2 = warm
recode warm2 1=1 2=3 3=2 4=4
slogit warm2 yr89 male white age ed prst, nolog
est store warm2
slogit warm yr89 male white age ed prst, nolog
est store warm
estimates table warm warm2, b(%9.3f) t(%6.2f)

* dimensions in slogit
slogit warm yr89 male white age ed prst, nolog
slogit warm yr89 male white age ed prst, nolog dim(2)

log close