Laboratorio 5
Federico Daverio
rm(list=ls())
setwd("C:/DAVE2/CIDE/Doc/Econometria II/Lab_5")
options(scipen=999)
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.0.5## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.4 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1## Warning: package 'ggplot2' was built under R version 4.0.5## Warning: package 'tibble' was built under R version 4.0.5## Warning: package 'tidyr' was built under R version 4.0.5## Warning: package 'readr' was built under R version 4.0.5## Warning: package 'dplyr' was built under R version 4.0.5## Warning: package 'forcats' was built under R version 4.0.5## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()## Warning: package 'janitor' was built under R version 4.0.5##
## Attaching package: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test## Warning: package 'modelsummary' was built under R version 4.0.5## Loading required package: dfidx## Warning: package 'dfidx' was built under R version 4.0.5##
## Attaching package: 'dfidx'## The following object is masked from 'package:stats':
##
## filter## Warning: package 'sandwich' was built under R version 4.0.5## Loading required package: msmModelos de conteo
Ejemplo 1. El número de personas muertas por patadas de mula o caballo en el ejército prusiano por año. Ladislaus Bortkiewicz recopiló datos de 20 volúmenes de Preussischen Statistik. Estos datos se recopilaron en 10 cuerpos del ejército prusiano a fines del siglo XIX a lo largo de 20 años.
Ejemplo 2. La cantidad de personas en fila frente a usted en la tienda de comestibles. Los predictores pueden incluir la cantidad de artículos que se ofrecen actualmente a un precio con descuento especial y si faltan tres días o menos para un evento especial (por ejemplo, un día festivo, un gran evento deportivo).
Ejemplo 3. El número de premios obtenidos por los estudiantes en una escuela secundaria. Los predictores de la cantidad de premios obtenidos incluyen el tipo de programa en el que se inscribió el estudiante (por ejemplo, vocacional, general o académico) y la puntuación en su examen final de matemáticas.
Descripciòn Datos
- n_award = numeros premios, variable dependiente
- math = variable continua de un puntaje en un test de matematica
- prog = programa de que hacen parte: 1 = “General”, 2 = “Academic” and 3 = “Vocational”
p <- read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")
p <- within(p, {
prog <- factor(prog, levels=1:3, labels=c("General", "Academic",
"Vocational"))
id <- factor(id)
})
summary(p)## id num_awards prog math
## 1 : 1 Min. :0.00 General : 45 Min. :33.00
## 2 : 1 1st Qu.:0.00 Academic :105 1st Qu.:45.00
## 3 : 1 Median :0.00 Vocational: 50 Median :52.00
## 4 : 1 Mean :0.63 Mean :52.65
## 5 : 1 3rd Qu.:1.00 3rd Qu.:59.00
## 6 : 1 Max. :6.00 Max. :75.00
## (Other):194Tenemos 200 observaciones. La media incondicional y la varianza no cambian mucho. Nuestras hipotesis dicen que estos valores condicionados con respecto a las observables relevantes sean similares.
Vamos ahora a anlizar estos dos momentos centrales de nuestras distribuciones en funciòn del programa de que hacen parte los estudaintes:
with(p, tapply(num_awards, prog, function(x) {
sprintf("M (SD) = %1.2f (%1.2f)", mean(x), sd(x))
}))## General Academic Vocational
## "M (SD) = 0.20 (0.40)" "M (SD) = 1.00 (1.28)" "M (SD) = 0.24 (0.52)"Y ploteamos un histograma para que resulte màs claro:
Que metodos de analisis deberìamos considerar?
Regresión de Poisson: la regresión de Poisson se utiliza a menudo para modelar datos de conteo Tiene extensiones utiles.
Regresión binomial negativa: la regresión binomial negativa se puede utilizar para datos de recuento excesivamente dispersos, es decir, cuando la varianza condicional excede la media condicional. Se puede considerar como una generalización de la regresión de Poisson ya que tiene la misma estructura media que la regresión de Poisson y tiene un parámetro adicional para modelar la sobredispersión. Si la distribución condicional de la variable de resultado está demasiado dispersa, es probable que los intervalos de confianza para la regresión binomial negativa sean más estrechos en comparación con los de una regresión de Poisson.
Modelo de regresión inflado con cero: los modelos inflados con cero intentan tener en cuenta el exceso de ceros. En otras palabras, se cree que existen dos tipos de ceros en los datos, “ceros verdaderos” y “ceros en exceso”. Los modelos inflados con cero estiman dos ecuaciones simultáneamente, una para el modelo de conteo y otra para el exceso de ceros.
Regresión OLS: las variables de resultado del recuento a veces se transforman logarítmicamente y se analizan mediante regresión OLS. Surgen muchos problemas con este enfoque, incluida la pérdida de datos debido a valores indefinidos generados al tomar el logaritmo de cero (que no está definido) y estimaciones sesgadas.
Utilizamos la funciòn glm para obtener nuestras desimaciones de un modelo de Poisson:
##
## Call:
## glm(formula = num_awards ~ prog + math, family = "poisson", data = p)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2043 -0.8436 -0.5106 0.2558 2.6796
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.24712 0.65845 -7.969 0.0000000000000016 ***
## progAcademic 1.08386 0.35825 3.025 0.00248 **
## progVocational 0.36981 0.44107 0.838 0.40179
## math 0.07015 0.01060 6.619 0.0000000000362501 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 287.67 on 199 degrees of freedom
## Residual deviance: 189.45 on 196 degrees of freedom
## AIC: 373.5
##
## Number of Fisher Scoring iterations: 6Cameron y Trivedi (2009) recomendaron el uso de errores estándar robustos para las estimaciones de los parámetros para controlar la violación leve del supuesto de distribución de que la varianza es igual a la media. Usamos el paquete R sandwich a continuación para obtener los errores estándar robustos.
cov.m1 <- vcovHC(m1, type="HC0")
std.err <- sqrt(diag(cov.m1))
r.est <- cbind(Estimate= coef(m1), "Robust SE" = std.err,
"Pr(>|z|)" = 2 * pnorm(abs(coef(m1)/std.err), lower.tail=FALSE),
LL = coef(m1) - 1.96 * std.err,
UL = coef(m1) + 1.96 * std.err)Analizando la salida del GLM noitamos que:
Los residuos de desviación se distribuyen aproximadamente normalmente si el modelo se especifica correctamente. En nuestro ejemplo, muestra un poco de asimetría ya que la mediana no es del todo cero.
El coeficiente de matemáticas es .07. Esto significa que el recuento de registros esperado para un aumento de una unidad en matemáticas es .07.
La variable indicadora progAcademic se compara entre prog = “Académico” y prog = “General”, el recuento de registros esperado para prog = “Académico” aumenta en aproximadamente 1,1. La variable indicadora prog.Vocational es la diferencia esperada en el recuento de registros ((aproximadamente .37)) entre prog = “Vocational” y el grupo de referencia (prog = “General”).
Podemos usar la desviación residual para realizar una prueba de bondad de ajuste para el modelo general. La desviación residual es la diferencia entre la desviación del modelo actual y la desviación máxima del modelo ideal donde los valores predichos son idénticos a los observados. Por lo tanto, si la diferencia residual es lo suficientemente pequeña, la prueba de bondad de ajuste no será significativa, lo que indica que el modelo se ajusta a los datos. Concluimos que el modelo se ajusta razonablemente bien porque la prueba de bondad de ajuste de chi-cuadrado no es estadísticamente significativa. Si la prueba hubiera sido estadísticamente significativa, indicaría que los datos no se ajustan bien al modelo.
with(m1, cbind(res.deviance = deviance, df = df.residual,
p = pchisq(deviance, df.residual, lower.tail=FALSE)))## res.deviance df p
## [1,] 189.4496 196 0.6182274También podemos probar el efecto general del prog comparando la desviación del modelo completo con la desviación del modelo excluyendo el prog. La prueba de chi-cuadrado de dos grados de libertad indica que el prog, en conjunto, es un predictor estadísticamente significativo de num_awards.
## Analysis of Deviance Table
##
## Model 1: num_awards ~ math
## Model 2: num_awards ~ prog + math
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 198 204.02
## 2 196 189.45 2 14.572 0.0006852 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A veces, es posible que deseemos presentar los resultados de la regresión como razones de tasas de incidencia y sus errores estándar, junto con el intervalo de confianza. Para calcular el error estándar de las tasas de incidencia de incidentes, usaremos el método Delta. Para ello, utilizamos la función deltamethod implementada en el paquete R msm.
s <- deltamethod(list(~ exp(x1), ~ exp(x2), ~ exp(x3), ~ exp(x4)),
coef(m1), cov.m1)
## exponentiate old estimates dropping the p values
rexp.est <- exp(r.est[, -3])
## replace SEs with estimates for exponentiated coefficients
rexp.est[, "Robust SE"] <- s
rexp.est## Estimate Robust SE LL UL
## (Intercept) 0.00526263 0.00339965 0.001483604 0.01866757
## progAcademic 2.95606545 0.94903937 1.575550534 5.54620292
## progVocational 1.44745846 0.57958742 0.660334536 3.17284023
## math 1.07267164 0.01119351 1.050955210 1.09483681El resultado anterior indica que la tasa de incidentes para prog = “Académico” es 2,96 veces la tasa de incidentes para el grupo de referencia (prog = “General”). Asimismo, la tasa de incidencia para prog = “Vocacional” es 1,45 veces la tasa de incidencia para el grupo de referencia que mantiene constantes las demás variables. El cambio porcentual en la tasa de incidentes de num_awards es del 7% por cada unidad de aumento en matemáticas.
Con el modelo de Poisson multiplicativo, los exponentes de los coeficientes son iguales a la razón de la tasa de incidencia (riesgo relativo).
Si queremos analizar los efectos marginales utilizaremos la funciòn predict.
Acomodamos los datos con un valor promedio de math score
## math prog
## 1 52.645 General
## 2 52.645 Academic
## 3 52.645 VocationalAplicamos el predict:
## $fit
## 1 2 3
## 0.2114109 0.6249446 0.3060086
##
## $se.fit
## 1 2 3
## 0.07050108 0.08628117 0.08833706
##
## $residual.scale
## [1] 1En el resultado anterior, vemos que el número predicho de eventos para el nivel 1 de prog es aproximadamente .21, manteniendo las matemáticas en su media. El número predicho de eventos para el nivel 2 de prog es mayor a .62, y el número predicho de eventos para el nivel 3 de prog es aproximadamente .31. Las proporciones de estos conteos predichos ((frac {.625} {. 211} = 2.96), (frac {.306} {. 211} = 1.45)) coinciden con lo que vimos al observar la TIR.
Podemos graficar el número predicho de eventos. El gráfico indica que se predice la mayor cantidad de premios para aquellos en el programa académico (prog = 2), especialmente si el estudiante tiene un puntaje alto en matemáticas. El número más bajo de premios previstos es para aquellos estudiantes en el programa general (prog = 1).
p$phat <- predict(m1, type="response")
## order by program and then by math
p <- p[with(p, order(prog, math)), ]
## create the plot
ggplot(p, aes(x = math, y = phat, colour = prog)) +
geom_point(aes(y = num_awards), alpha=.5, position=position_jitter(h=.2)) +
geom_line(size = 1) +
labs(x = "Math Score", y = "Expected number of awards")
## Importante
Cuando parece haber un problema de dispersión, primero debemos verificar si nuestro modelo está especificado adecuadamente, controlando por variables omitidas y las formas funcionales. Por ejemplo, si omitimos la variable predictora prog en el ejemplo anterior, nuestro modelo parecería tener un problema de dispersión excesiva. En otras palabras, un modelo mal especificado podría presentar un síntoma como un problema de dispersión excesiva.
Suponiendo que el modelo se especifica correctamente, debe comprobarse el supuesto de que la varianza condicional es igual a la media condicional. Hay varias pruebas que incluyen la prueba de razón de verosimilitud del parámetro alfa de sobredispersión ejecutando el mismo modelo utilizando una distribución binomial negativa. El paquete R pscl (Laboratorio Computacional de Ciencias Políticas, Universidad de Stanford) proporciona muchas funciones para datos binomiales y de conteo, incluido odTest para probar la sobredispersión.
Una causa común de dispersión excesiva es el exceso de ceros, que a su vez se generan mediante un proceso de generación de datos adicional. En esta situación, se debe considerar el modelo inflado cero.
Si el proceso de generación de datos no permite ningún 0 (como el número de días en el hospital), entonces un modelo truncado en cero puede ser más apropiado.
La variable de resultado en una regresión de Poisson no puede tener números negativos y la exposición no puede tener ceros.
Existen muchas medidas diferentes de pseudo-R-cuadrado. Todos intentan proporcionar información similar a la proporcionada por R-cuadrado en la regresión MCO, aunque ninguno de ellos puede interpretarse exactamente como se interpreta R-cuadrado en la regresión MCO.
Regresiòn Binomia Negativa
## Loading required package: foreign## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dfidx':
##
## select## The following object is masked from 'package:dplyr':
##
## selectEjemplos
Ejemplo 1. Los administradores escolares estudian el comportamiento de asistencia de los estudiantes de tercer año de secundaria en dos escuelas. Los predictores del número de días de ausencia incluyen el tipo de programa en el que está matriculado el estudiante y una prueba estandarizada de matemáticas.
Ejemplo 2. Un investigador relacionado con la salud está estudiando el número de visitas al hospital en los últimos 12 meses por parte de personas mayores en una comunidad en función de las características de los individuos y los tipos de planes de salud bajo los cuales cada uno está cubierto.
Datos
Tenemos datos de asistencia de 314 estudiantes de tercer año de secundaria de dos escuelas secundarias urbanas en el archivo nb_data. La variable de respuesta de interés es días ausentes, días ab. La variable matemática proporciona el puntaje matemático estandarizado para cada estudiante. La variable prog es una variable nominal de tres niveles que indica el tipo de programa educativo en el que está matriculado el estudiante.
dat <- read.dta("https://stats.idre.ucla.edu/stat/stata/dae/nb_data.dta")
dat <- within(dat, {
prog <- factor(prog, levels = 1:3, labels = c("General", "Academic", "Vocational"))
id <- factor(id)
})
summary(dat)## id gender math daysabs prog
## 1001 : 1 female:160 Min. : 1.00 Min. : 0.000 General : 40
## 1002 : 1 male :154 1st Qu.:28.00 1st Qu.: 1.000 Academic :167
## 1003 : 1 Median :48.00 Median : 4.000 Vocational:107
## 1004 : 1 Mean :48.27 Mean : 5.955
## 1005 : 1 3rd Qu.:70.00 3rd Qu.: 8.000
## 1006 : 1 Max. :99.00 Max. :35.000
## (Other):308Hacemos tambièn un analisis graficas:
ggplot(dat, aes(daysabs, fill = prog)) + geom_histogram(binwidth = 1) + facet_grid(prog ~
., margins = TRUE, scales = "free")
Cada variable tiene 314 observaciones válidas y sus distribuciones parecen bastante razonables. La media incondicional de nuestra variable de resultado es mucho menor que su varianza.
La siguiente tabla muestra el número promedio de días ausentes por tipo de programa y parece sugerir que el tipo de programa es un buen candidato para predecir el número de días ausentes, nuestra variable de resultado, porque el valor medio del resultado parece variar por prog. Las variaciones dentro de cada nivel de prog son más altas que las medias dentro de cada nivel. Estos son los promedios y variaciones condicionales. Estas diferencias sugieren que existe una dispersión excesiva y que sería apropiado un modelo binomial negativo.
## General Academic Vocational
## "M (SD) = 10.65 (8.20)" "M (SD) = 6.93 (7.45)" "M (SD) = 2.67 (3.73)"Nota a margen
La función tapply aplica (de ahí parte de su nombre) una función a un vector en los subvectores que define otro vector máscara.
tapply(iris$Petal.Length, iris$Species, mean)
En el ejemplo anterior, tapply aplica la función mean a un vector, la longitud del pétalo, pero de especie en especie. Responde por tanto a la pregunta de cuál es la media de la longitud del pétalo por especie. Es una operación similar a la conocida en SQL como group by
tapply es una de las llamadas funciones de orden superior porque acepta como argumento otra función (mean en el caso anterior). Las funciones a las que llaman tanto tapply como el resto de las funciones de orden superior pueden tener parámetros adicionales. Por ejemplo, si el vector sobre el que se quiere calcular la media contiene NA, el resultado es NA.
Análisis de regresión binomial negativa
A continuación, usamos la función glm.nb del paquete MASS para estimar una regresión binomial negativa.
##
## Call:
## glm.nb(formula = daysabs ~ math + prog, data = dat, init.theta = 1.032713156,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1547 -1.0192 -0.3694 0.2285 2.5273
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.615265 0.197460 13.245 < 0.0000000000000002 ***
## math -0.005993 0.002505 -2.392 0.0167 *
## progAcademic -0.440760 0.182610 -2.414 0.0158 *
## progVocational -1.278651 0.200720 -6.370 0.000000000189 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.0327) family taken to be 1)
##
## Null deviance: 427.54 on 313 degrees of freedom
## Residual deviance: 358.52 on 310 degrees of freedom
## AIC: 1741.3
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.033
## Std. Err.: 0.106
##
## 2 x log-likelihood: -1731.258R primero muestra la llamada y los residuales de desviación. A continuación, vemos los coeficientes de regresión para cada una de las variables, junto con los errores estándar, las puntuaciones z y los valores p. La variable matemática tiene un coeficiente de -0,006, que es estadísticamente significativo. Esto significa que por cada aumento de una unidad en matemáticas, el recuento de registros esperado del número de días ausentes disminuye en 0,006. La variable indicadora que se muestra como progAcademic es la diferencia esperada en el recuento de registros entre el grupo 2 y el grupo de referencia (prog = 1). El recuento de registros esperado para el nivel 2 de prog es 0.44 menor que el recuento de registros esperado para el nivel 1. La variable indicadora de progVocational es la diferencia esperada en el recuento de registros entre el grupo 3 y el grupo de referencia. es 1,28 menor que el recuento de registros esperado para el nivel 1.
Para determinar si el prog en sí, en general, es estadísticamente significativo, podemos comparar un modelo con y sin prog. La razón por la que es importante ajustar modelos separados es que, a menos que lo hagamos, el parámetro de sobredispersión se mantiene constante.
## Likelihood ratio tests of Negative Binomial Models
##
## Response: daysabs
## Model theta Resid. df 2 x log-lik. Test df LR stat.
## 1 math 0.8558565 312 -1776.306
## 2 math + prog 1.0327132 310 -1731.258 1 vs 2 2 45.04798
## Pr(Chi)
## 1
## 2 0.000000000165179La prueba de chi-cuadrado de dos grados de libertad indica que el prog es un predictor estadísticamente significativo de daysabs.
El parámetro theta que se muestra es el parámetro de dispersión. En R pero es el inverso del parametro original. Con cuidado.
Averiguar hipotesis
Los modelos binomiales negativos asumen que las medias condicionales no son iguales a las varianzas condicionales. Esta desigualdad se captura estimando un parámetro de dispersión (que no se muestra en el resultado) que se mantiene constante en un modelo de Poisson. Por lo tanto, el modelo de Poisson en realidad está anidado en el modelo binomial negativo. Luego, podemos usar una prueba de razón de verosimilitud para comparar estos dos y probar esta suposición del modelo. Para hacer esto, ejecutaremos nuestro modelo como Poisson.
m3 <- glm(daysabs ~ math + prog, family = "poisson", data = dat)
pchisq(2 * (logLik(m1) - logLik(m3)), df = 1, lower.tail = FALSE)## 'log Lik.' 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002157298 (df=5)En este ejemplo, el valor asociado de chi-cuadrado estimado es 926,03 con un grado de libertad. Esto sugiere fuertemente que el modelo binomial negativo, que estima el parámetro de dispersión, es más apropiado que el modelo de Poisson.
Podemos obtener los intervalos de confianza para los coeficientes perfilando la función de verosimilitud.
## Waiting for profiling to be done...## Estimate 2.5 % 97.5 %
## (Intercept) 2.615265446 2.24205576 3.012935926
## math -0.005992988 -0.01090086 -0.001066615
## progAcademic -0.440760012 -0.81006586 -0.092643481
## progVocational -1.278650721 -1.68348970 -0.890077623Podríamos estar interesados en observar las tasas de incidencia en lugar de los coeficientes. Para hacer esto, podemos exponencializar los coeficientes de nuestro modelo. Lo mismo se aplica a los intervalos de confianza.
## Estimate 2.5 % 97.5 %
## (Intercept) 13.6708448 9.4126616 20.3470498
## math 0.9940249 0.9891583 0.9989340
## progAcademic 0.6435471 0.4448288 0.9115184
## progVocational 0.2784127 0.1857247 0.4106239El resultado anterior indica que la tasa de incidentes para prog = 2 es 0,64 veces la tasa de incidentes para el grupo de referencia (prog = 1). Asimismo, la tasa de incidentes para prog = 3 es 0,28 veces la tasa de incidentes para el grupo de referencia que mantiene constantes las demás variables. El cambio porcentual en la tasa de incidencia de daysabs es una disminución del 1% por cada unidad de aumento en matemáticas.
La forma de la ecuación del modelo para la regresión binomial negativa es la misma que la de la regresión de Poisson. El registro del resultado esperado se predice con una combinación lineal de predictores:
\[ ln(\widehat{daysabs_i}) = Intercept + b_1I(prog_i = 2) + b_2I(prog_i = 3) + b_3math_i \]
Por lo tanto tendremos que:
\[\widehat{daysabs_i} = e^{Intercept + b_1 I(prog_i = 2) + b_2I(prog_i = 3) + b_3 math_i} = e^{Intercept}e^{b_1 I(prog_i = 2)}e^{b_2 I(prog_i = 3)}e^{b_3 math_i}\]
Entocnes los coeficientes tienen un efecto aditivo en la escala \(ln(y)\) y la TIR tiene un efecto multiplicativo en la escala \(y\). El parámetro de dispersión en la regresión binomial negativa no afecta los recuentos esperados, pero sí afecta la varianza estimada de los recuentos esperados.
Valores predecidos
Para calcular efectos marginales podemos:
newdata1 <- data.frame(math = mean(dat$math), prog = factor(1:3, levels = 1:3,
labels = levels(dat$prog)))
newdata1$phat <- predict(m1, newdata1, type = "response")
newdata1## math prog phat
## 1 48.26752 General 10.236899
## 2 48.26752 Academic 6.587927
## 3 48.26752 Vocational 2.850083En el resultado anterior, vemos que el número predicho de eventos (por ejemplo, días ausentes) para un programa general es de aproximadamente 10,24, manteniendo las matemáticas en su media. El número previsto de eventos para un programa académico es inferior a 6,59 y el número previsto de eventos para un programa vocacional es de aproximadamente 2,85.
A continuación, obtendremos el número medio predicho de eventos para valores de matemáticas en todo su rango para cada nivel de prog y graficaremos estos.
newdata2 <- data.frame(
math = rep(seq(from = min(dat$math), to = max(dat$math), length.out = 100), 3),
prog = factor(rep(1:3, each = 100), levels = 1:3, labels =
levels(dat$prog)))
newdata2 <- cbind(newdata2, predict(m1, newdata2, type = "link", se.fit=TRUE))
newdata2 <- within(newdata2, {
DaysAbsent <- exp(fit)
LL <- exp(fit - 1.96 * se.fit)
UL <- exp(fit + 1.96 * se.fit)
})
ggplot(newdata2, aes(math, DaysAbsent)) +
geom_ribbon(aes(ymin = LL, ymax = UL, fill = prog), alpha = .25) +
geom_line(aes(colour = prog), size = 2) +
labs(x = "Math Score", y = "Predicted Days Absent")
El gráfico muestra el recuento esperado en el rango de puntajes matemáticos, para cada tipo de programa junto con intervalos de confianza del 95 por ciento. Tenga en cuenta que las líneas no son rectas porque este es un modelo logarítmico lineal, y lo que se traza son los valores esperados, no el logaritmo de los valores esperados.
Cuidados
- No se recomienda aplicar modelos binomiales negativos a muestras pequeñas.
- Una causa común de dispersión excesiva es el exceso de ceros por un proceso de generación de datos adicional. En esta situación, se debe considerar el modelo inflado cero.
- Si el proceso de generación de datos no permite ningún 0 (como el número de días en el hospital), entonces un modelo truncado en cero puede ser más apropiado.
#Referencias: 1. https://stats.idre.ucla.edu/r/dae/poisson-regression/