log using ch4binary_1nov07, text replace // Regression Models for Categorical Dependent Variables - 2nd Edition // Chapter 4 - Models for Binary Outcomes LOGIT/PROBIT // Long and Freese - 27Jul2005 // Comentarios en español: Javier Aparicio - 27oct2008 *** CONTENIDO *** Pruebas de hipótesis con MLE: test y lrtest *** Análisis de residuales y detección de valores atípicos (outliers): standardized residuals, *** Cook's distance, lowess. *** Medidas de bondad de ajuste: fitstat *** Analizando las probabilidades predichas por LOGIT *** Comandos de SPOST: fitstat, leastlikely, prvalue, prtab, prgen, prchange, listcoef /* Para ejecutar este do-file necesitas instalar las rutinas de SPOST así como instalar las bases de datos de los ejemplos. En stata9, usar: findit spost9_ado // instala los comandos y files de Spost findit spost9_do // copia al dir de trabajo bases de datos y do-files con ejemplos del libro */ * Algunas opciones iniciales version 9 // fija la versión de stata9 (en caso de usar stata10) *set scheme s2manual // fija las gráficas en estilo monocromático *set more off // para no hacer pausas en cada pantalla de salida *** LOGIT y PROBIT // estimation using logit and probit (pg 115 RevEd) sysuse binlfp2, clear //si el comando sysuse no funciona, intenta "use..." desc lfp k5 k618 age wc hc lwg inc summarize lfp k5 k618 age wc hc lwg inc logit lfp k5 k618 age wc hc lwg inc, nolog estimates store logit probit lfp k5 k618 age wc hc lwg inc, nolog estimates store probit * Comparing the estimates from logit and probit estimates table logit probit, b(%9.3f) t label varwidth(30) *** Pruebas de hipótesis con MLE logit lfp k5 k618 age wc hc lwg inc, nolog * 1. Wald test test k5 display sqrt(55.14) display sqrt(r(chi2)) // or, use returns to get the value of 55.14 * 2. Likelihood ratio test (LR test) - es similar a la prueba F de OLS logit lfp k5 k618 age wc hc lwg inc, nolog estimates store fmodel // guardamos los resultados del "modelo completo" logit lfp k618 age wc hc lwg inc, nolog estimates store nmodel // guardamos resultados del modelo "restringido o anidado" lrtest fmodel nmodel // comparamos ambos lrtest fmodel . // el . se refiere al último modelo en memoria // Testing multiple coefficients (pg 122 RevEd) logit lfp k5 k618 age wc hc lwg inc, nolog * Wald tests test hc wc // Ho: beta_hc y beta_wc son "conjuntamente insignificantes" test hc = wc // Ho: beta_hc es igual a beta_wc * Likelihood ratio tests * LR tests - similares a la prueba F de significancia conjunta en OLS logit lfp k5 k618 age wc hc lwg inc, nolog // el modelo "completo" o no restringido estimates store fullmodel logit lfp k5 k618 age lwg inc, nolog // una version restringida o "anidada" del fullmodel estimates store nestedmodel lrtest fullmodel nestedmodel // Ho: la diferencia en log-likelihood entre ambos modelos // es estadísticamente insignificante logit lfp, nolog // el modelo más restringido posible: solo con intercepto estimates store intercept_only lrtest fullmodel intercept_only * esta comparación entre el modelo completo con el más restringido es el LR test * que aparece en el output estándar de logit *** Análisis de residuales y detección de valores atípicos (outliers) // Analizando los residuales de logit (pg 125 RevEd) logit lfp k5 k618 age wc hc lwg inc, nolog predict rstd, rst // obtengo los residuales estandarizados del modelo label var rstd "Pearson Standardized Residual" sort inc, stable // ordenamos los datos de acuerdo a ingreso generate index = _n // generamos un "índice" de 1 hasta n label var index "Observation Number" * index plot of std pearson residuals graph twoway scatter rstd index, xlabel(0(200)800) ylabel(-4(2)4) /// xtitle("Observation Number") yline(0) msymbol(Oh) /// ysize(2.7051) xsize(4.0413) * residual = Yobs - Ypredicha... * si residual > 0 -> Yobs=1 mientras que Ypredicha=0, y viceversa: * si residual < 0 -> Yobs=0 mientras que Ypredicha=1 graph export 04rstd.emf , replace * guarda la gráfica en el dir. de trabajo (hay varios formatos, ver ayuda) * En las gráficas puedes usar etiquetas o valores en vez de puntitos: * index plot of std pearson residuals labeled with the index graph twoway scatter rstd index, xlabel(0(200)800) ylabel(-4(2)4) /// xtitle("Observation Number") yline(0) /// msymbol(none) mlabel(index) mlabposition(0) /// ysize(2.7051) xsize(4.0413) graph export 04rstdcase.emf , replace * index plot of std pearson residuals labeled with value of k5 graph twoway scatter rstd index, /// msymbol(none) mlabel(k5) mlabposition(0) /// caption("Values indicate # of young children") /// xlabel(0(200)800) xtitle("Observation Number") /// ylabel(-4(2)4) yline(0) ysize(2.7051) xsize(4.0413) * Graficando solamente las observacione más atípicas (ie, con residuales grandes) * index plot of std pearson residuals for large residuals only (atypical values) * ie, std residuals > |1.7| std. dev. graph twoway scatter rstd index if (rstd>1.7) | (rstd<-1.7), /// msymbol(none) mlabel(k5) mlabposition(0) /// caption("Values indicate # of young children") /// xlabel(0(200)800) xtitle("Observation Number") /// ylabel(-4(2)4) yline(0) ysize(2.7051) xsize(4.0413) graph export 04rstdbig.emf , replace * list a single point list in 142, noobs // la observacion #142 ordenada según ingreso * list large outliers (std residual > | 2.5 std. dev.) list rstd index if rstd>2.5 | rstd<-2.5 // Identifying influential cases (pg 128 RevEd) predict cook, dbeta // see "help logit postestimation" label var cook "Cook's Statistic" * cook/dbeta calculates the "Pregibon Delta-Beta" influence statistic, a standardized measure of * the difference in the coefficient vector due to deletion of the observation along * with all others that share the same covariate pattern. graph twoway scatter cook index, xlabel(0(200)800) ylabel(0(.1).3) /// xtitle("Observation Number") yline(.1 .2) /// msymbol(none) mlabel(index) mlabposition(0) /// ysize(2.7051) xsize(4.0413) graph export 04cookcase.emf , replace // Identifying "least likely" cases // Cuáles son los casos observados como 1's (o 0's) que tienen la menor probabilidad predicha // de haberlo sido según el modelo estimado. logit lfp k5 k618 age wc hc lwg inc, nolog leastlikely k5 k618 wc inc * Por ejemplo, el caso #752 tiene una muy baja probabilidad predicha de trabajar (9%), dado que * tiene cuatro hijos y un muy elevado ingreso familiar... y sin embargo trabaja. *** Bondad de ajuste (goodness of fit) *** en MLE no existen residuales cuadrados ni R2, así que hay una variadad de fit statistics // scalar measures of fit (pg 129 RevEd) logit lfp k5 k618 age wc hc lwg inc, nolog estimates store model1 // guardamos los estimadores de model1 fitstat, saving(m1) // muestra y guarda diferentes estadísticos de ajuste logit lfp age wc hc lwg inc, nolog // modelo1 pero sin kids estimates store model2 estimates table model1 model2, b(%9.3f) t fitstat, using(m1) // compara el modelo en memoria con m1 * alternativamente quietly logit lfp k5 k618 age wc hc lwg inc, nolog quietly fitstat, save quietly logit lfp age wc hc lwg inc, nolog fitstat, diff // compara modelo en memoria con modelo guardado // hosmer-lemeshow test (pg 2Ed) logit lfp k5 k618 age wc hc lwg inc, nolog estat gof, group(10) table estat gof * lowess plot predict p1 // guarda las probabilidades predichas del modelo en la variable p1 lowess lfp p1, ylabel(0(.2)1, grid) xlabel(0(.2)1, /// grid) ysize(5) xsize(5) addplot(function y = x, legend(off)) //addplot añade la recta de 45º * lowess calcula y grafica una "locally weighted regression" de Yvar sobre Xvar y la grafica. * en este caso, dada la frecuencia relativa de 1's y 0's (en el eje Y), se calcula una * "proporción local" que se compara con la probabilidad predicha por el modelo (eje X) graph export 04hllowess.emf, replace * Si la curva lowess está muy cerca de la recta de 45º, el modelo tiene "buen ajuste", y * allí donde la curva lowess se aleja de la recta de 45º, el modelo tiene un peor ajuste. * En este caso, el modelo tiene un pobre ajuste cuando p1 es bajo (ie, bajas probabilidades * de trabajar). *** Analizando las probabilidades predichas por LOGIT // predicted probabilities with predict (pg 132 RevEd) * predictions from logit sysuse binlfp2, clear logit lfp k5 k618 age wc hc lwg inc, nolog predict prlogit summarize prlogit label var prlogit "Logit: Pr(lfp)" dotplot prlogit, ylabel(0(.2)1) /// ysize(2.7051) xsize(4.0413) graph export 04dotpredict.emf, replace * comparing logit and probit predictions sysuse binlfp2, clear logit lfp k5 k618 age wc hc lwg inc, nolog predict prlogit label var prlogit "Logit: Pr(lfp)" probit lfp k5 k618 age wc hc lwg inc, nolog predict prprobit label var prprobit "Probit: Pr(lfp)" pwcorr prlogit prprobit * graphing predicted probabilities from logit and probit graph twoway scatter prlogit prprobit, /// xlabel(0(.25)1) ylabel(0(.25)1) /// xline(.25(.25)1) yline(.25(.25)1) /// plotregion(margin(zero)) msymbol(Oh) /// ysize(4.0413) xsize(4.0413) graph export 04logitprobit.emf, replace ** PRVALUE - predicción de probabilidad de casos particulares // individual predicted probabilities with prvalue (pg 135 RevEd) logit lfp k5 k618 age wc hc lwg inc, nolog * young, low income, low education families with young children. prvalue, x(age=35 k5=2 wc=0 hc=0 inc=15) rest(mean) * highly educated families with no children at home. prvalue, x(age=50 k5=0 k618=0 wc=1 hc=1) rest(mean) * while an average person... prvalue, rest(mean) ** PRTAB - tablas de probabilidad predicha // tables of predicted probabilities with prtab (pg 136 RevEd) logit lfp k5 k618 age wc hc lwg inc, nolog prtab k5 wc, rest(mean) // ojo, sólo tabula variables discretas * we can use prvalue to get the same information plus confidence intervals prvalue, x(k5=0 wc=0) rest(mean) brief prvalue, x(k5=1 wc=0) rest(mean) brief prvalue, x(k5=2 wc=0) rest(mean) brief prvalue, x(k5=3 wc=0) rest(mean) brief * compute the differences: disp 0.6069-0.7758 disp 0.2633-0.4449 disp 0.0764-0.1565 disp 0.0188-0.0412 ** PRGEN - genera "series" de probabilidades predichas para hacer graficas // graphing predicted probabilities with -prgen- (pg 137 RevEd) * compute predictions at age 30 prgen inc, from(0) to(100) generate(p30) x(age=30) rest(mean) ncases(11) label var p30p1 "Age 30" label var p30p0 "Age 30" * Calcula la prob. predicha conforme family income va de 0 a 100, con 11 intervalos * (0, 10, 20 ... 100), manteniendo edad=30 y el resto de variables en sus medias * ie, individuos de 30 años con características promedio y diferentes niveles de ingreso * Noten que se generan 3 variables: p30x (la escala de edad) p30p0 (prob Y=0) y p30p1 (prob Y=1) * compute predictions at age 40 prgen inc, from(0) to(100) generate(p40) x(age=40) rest(mean) n(11) label var p40p1 "Age 40" label var p40p0 "Age 40" * compute predictions at age 50 prgen inc, from(0) to(100) generate(p50) x(age=50) rest(mean) n(11) label var p50p1 "Age 50" label var p50p0 "Age 50" * compute predictions at age 60 prgen inc, from(0) to(100) generate(p60) x(age=60) rest(mean) n(11) label var p60p1 "Age 60" label var p60p0 "Age 60" * list and graphing the predictions list p30p1 p40p1 p50p1 p60p1 p60x in 1/11 * Gráfica para p1 - probabilidad de SI trabajar graph twoway connected p30p1 p40p1 p50p1 p60p1 p60x, /// ytitle("Pr(In Labor Force)") ylabel(0(.25)1) /// xtitle("Income") /// ysize(2.7051) xsize(4.0413) graph export 04ageincome.emf, replace * Gráfica para p0 - probabilidad de NO trabajar graph twoway connected p30p0 p40p0 p50p0 p60p0 p60x, /// ytitle("Pr(NOT in Labor Force)") ylabel(0(.25)1) /// xtitle("Income") /// ysize(2.7051) xsize(4.0413) ** Graficando intervalos de confianza para probabilidades predichas // plotting confidence intervals prgen age, from(20) to(70) generate(prlfp) rest(mean) gap(2) ci * Calcula la prob. predicha para edades de 20 a 70 con intervalos de 2 años, * manteniendo todo lo demás en sus medias, y con intervalos de confianza. * Noten que se generan 7 variables: * prlfpx (escala de edad), prlfpp0 prlfpp1 (prob0 y prob1, respectivamente), * prlfpp0lb prlfpp1lb (95% lower bounds para p0 y p1) * prlfpp0ub prlfpp1up (95% upper bounds para p0 y p1) label var prlfpp1 "Predicted probability" label var prlfpp1ub "95% upper limit" label var prlfpp1lb "95% lower limit" label var prlfpx "Age" * using lines to show confidence interval twoway /// (connected prlfpp1 prlfpx, /// clcolor(black) clpat(solid) clwidth(medthick) /// msymbol(i) mcolor(none)) /// (connected prlfpp1ub prlfpx, /// msymbol(i) mcolor(none) /// clcolor(black) clpat(dash) clwidth(thin)) /// (connected prlfpp1lb prlfpx, /// msymbol(i) mcolor(none) /// clcolor(black) clpat(dash) clwidth(thin)), /// ytitle("Probability of Being in Labor Force") yscale(range(0 .35)) /// ylabel(, grid glwidth(medium) glpattern(solid)) /// xscale(range(20 70)) xlabel(20(10)70) /// ysize(2.7051) xsize(4.0413) graph export 04ageci.emf, replace * using shading to show confidence interval graph twoway /// (rarea prlfpp1lb prlfpp1ub prlfpx, color(gs14)) /// (connected prlfpp1 prlfpx, /// clcolor(black) clpat(solid) clwidth(medthick) /// msymbol(i) mcolor(none)), /// ytitle("Probability of Being in Labor Force") yscale(range(0 .35)) /// ylabel(, grid glwidth(medium) glpattern(solid)) /// xscale(range(20 70)) xlabel(20(10)70) /// ysize(2.7051) xsize(4.0413) /// legend(label (1 "95% confidence interval")) graph export 04ageci2.emf, replace ** PRGEN - Comparando probabilidades predichas para 2 o más grupos // prgen to compare those who do and do not attend college (not in book) prgen age, from(30) to(60) generate(wc1) x(wc=1) rest(mean) n(13) label var wc1p1 "Attended College" prgen age, from(30) to(60) generate(wc0) x(wc=0) rest(mean) n(13) label var wc0p1 "Did Not Attend College" graph twoway connected wc1p1 wc0p1 wc1x, /// xtitle("Age") ytitle("Pr(In Labor Force)") /// ylabel(0(.25)1) xlabel(30(10)60) ** PRCHANGE - calcular "cambios" en probabilidades predichas // changes in predicted probabilities (pg 139 RevEd) mfx compute, at(wc=1 age=40) // stata calcula efectos marginales * PRCHANGE cacula MAS efectos o cambios para una variable de interés: prchange age, x(wc=1 age=40) help // help añade notas de intrerpretación prchange, help // ...o para todas las variables prchange k5 age wc lwg inc, fromto // ...o para A, B y su diferencia * prchange no brinda intervalos de confianza, pero prvalue sí: prchange wc, brief // brief omite algunos datos de referencia prvalue, x(wc=0) rest(mean) save prvalue, x(wc=1) diff // compara con el prvalue salvado antes mfx compute // mfx tambien calcula el IC de wc * Calculando el prvalue para un grupo de variables binarias (en un loop) * automating prvalue for discrete change for 0/1 changes foreach v in k5 k618 wc { // inicia un ciclo di _n "** Change from 0 to 1 in `v'" qui prvalue, x(`v'=0) rest(mean) save prvalue, x(`v'=1) dif brief } // termina el ciclo * Calculando prvalue para cambios de 1 desv. est."centrada" * automating prvalue for centered sd changes foreach v in age lwg inc { qui sum `v' local start = r(mean) - (.5*r(sd)) local end = r(mean) + (.5*r(sd)) di _n "** Change from `start' to `end' in `v'" qui prvalue, x(`v'=`start') rest(mean) save prvalue, x(`v'=`end') dif brief } * discrete change using prvalue prvalue, x(age=30) save brief prvalue, x(age=40) diff brief * Cambios en prob. predicha para un cambio arbitrario en X (delta) * o bien un cambio arbitrario "no centrado" * discrete change using prchange with delta and uncentered options prchange age, x(age=30) delta(10) rest(mean) brief * o bien un cambio arbitrario "no centrado" prchange age, x(age=30) uncentered delta(10) rest(mean) brief ** LISTCOEF - lista coeficientes estimados en varios formatos // odds ratios using -listcoef- (pg 146 RevEd) quietly logit lfp k5 k618 age wc hc lwg inc, nolog listcoef, help * cuando beta es negativo, odds-ratio es menor a 1 (y viceversa) listcoef, reverse help // odds ratio invertidos listcoef, percent help // cambio porcentual en odds ratio listcoef, std help // betas estandarizadas log close