Recuerda que para enlaces canónicos, tenemos la igualdad entre el parámetro de la distribución y el componente sistemático, de modo que \(\theta_i=\eta _i=\mathbf{x}_i^{\mathbf{\prime }}\boldsymbol \beta\). Así, con \(\phi _i=\phi /w_i\), la función de log-verosimilitud es:
Ten en cuenta que, al igual que en las ecuaciones normales de la regresión lineal ordinaria, no necesitamos considerar la estimación del parámetro de escala de la varianza \(\phi\) en esta etapa. Es decir, primero podemos calcular \(\mathbf{b}_{MLE}\) y, cuando sea necesario, estimar \(\phi\). (Para ciertas distribuciones como la binomial y la de Poisson, \(\phi\) es conocido y no requiere estimación).
Como se describe en la Sección 11.9, los estimadores de máxima verosimilitud son consistentes y tienen distribuciones normales para muestras grandes bajo condiciones generales. La inferencia de máxima verosimilitud proporciona un mecanismo para calcular esta distribución. A partir de las ecuaciones (13.6) y del suplemento técnico sobre la verosimilitud (Sección 11.9, ecuación 11.14), la matriz de información correspondiente es:
La inversa de la matriz de información es la matriz de varianza-covarianza para muestras grandes de \(\mathbf{b}_{MLE}\). Específicamente, la raíz cuadrada del \((j+1)\)-ésimo elemento diagonal de la inversa de esta matriz proporciona el error estándar para \(b_{j,MLE}\), el cual denotamos como \(se(b_{j,MLE})\). Las extensiones a enlaces generales son similares.
# AutoDat <- read.csv("CSVData/AutoCollision.csv", quote="", header=TRUE)
AutoDat <- read.table("CSVData/AutoCollision1.txt", header=TRUE)
# MODELO CON DOS FACTORES
model1 <- glm(Severity~factor(Age)+factor(Vehicle_Use),
control = glm.control(maxit = 50), weights=Claim_Count,
family=Gamma(link="log"), data = AutoDat)
#summary(model1)
#model1$coefficients
Agecoeff <- model1$coefficients[1:8]
Agecoeff[1]<- 0
EAgecoeff <- exp(Agecoeff)
#EAgecoeff
Vehcoeff <- model1$coefficients[9:11];
Vehcoeff <- c(Vehcoeff,0)
EVehcoeff <- exp(Vehcoeff)
#EVehcoeff
Premiums <- exp(model1$coefficients[1])*EAgecoeff %*%t(EVehcoeff)
Premiums1 <- round(Premiums, digits = 2)
AgeGroup <- dimnames( table(AutoDat$Age) )
AgeGroup1 <- AgeGroup
# AgeGroup1 <- gsub("\226", "-", AgeGroup1[[1]] , fixed = TRUE)
tableout <- cbind(AgeGroup1, Premiums1)
colnames(tableout) <- c("Grupo de Edad", "ConducirCorto", "ConducirLargo", "Placer", "Negocios")
TableGen1(TableData=tableout,
TextTitle='Estimación de la Severidad Esperada para un Plan Tarifario Multiplicativo',
Align='c', ColumnSpec=1:3, Digits = 2, BorderRight = 1,
ColWidth = ColWidth4 )
Tabla 13.1: Estimación de la Severidad Esperada para un Plan Tarifario Multiplicativo
|
Grupo de Edad
|
ConducirCorto
|
ConducirLargo
|
Placer
|
Negocios
|
|
c(“17?20”, “21?24”, “25?29”, “30?34”, “35?39”, “40?49”, “50?59”, “60+”)
|
322.17
|
265.56
|
254.9
|
419.07
|
|
c(“17?20”, “21?24”, “25?29”, “30?34”, “35?39”, “40?49”, “50?59”, “60+”)
|
320.66
|
264.31
|
253.7
|
417.1
|
|
c(“17?20”, “21?24”, “25?29”, “30?34”, “35?39”, “40?49”, “50?59”, “60+”)
|
297.26
|
245.02
|
235.19
|
386.66
|
|
c(“17?20”, “21?24”, “25?29”, “30?34”, “35?39”, “40?49”, “50?59”, “60+”)
|
284.85
|
234.8
|
225.37
|
370.53
|
|
c(“17?20”, “21?24”, “25?29”, “30?34”, “35?39”, “40?49”, “50?59”, “60+”)
|
229.37
|
189.06
|
181.47
|
298.35
|
|
c(“17?20”, “21?24”, “25?29”, “30?34”, “35?39”, “40?49”, “50?59”, “60+”)
|
248.15
|
204.54
|
196.33
|
322.78
|
|
c(“17?20”, “21?24”, “25?29”, “30?34”, “35?39”, “40?49”, “50?59”, “60+”)
|
251.95
|
207.67
|
199.34
|
327.72
|
|
c(“17?20”, “21?24”, “25?29”, “30?34”, “35?39”, “40?49”, “50?59”, “60+”)
|
246.47
|
203.16
|
195
|
320.6
|
Sobredispersión
Para algunos miembros de la familia exponencial lineal, como las distribuciones Bernoulli y Poisson, la varianza está determinada por la media. En contraste, la distribución normal tiene un parámetro de escala separado. Al ajustar modelos a datos con variables dependientes binarias o de conteo, es común observar que la varianza excede la anticipada por el ajuste de los parámetros de la media. Este fenómeno se conoce como sobredispersión. Existen varios modelos probabilísticos alternativos disponibles para explicar este fenómeno, dependiendo de la aplicación en cuestión. Ver la Sección 12.3 para un ejemplo y McCullagh y Nelder (1989) para un inventario más detallado.
Aunque llegar a un modelo probabilístico satisfactorio es el camino más deseable, en muchas situaciones los analistas se conforman con postular un modelo aproximado mediante la relación:
\[
\mathrm{Var~}y_i=\sigma ^{2}\phi ~b^{\prime \prime }(\theta_i)/w_i.
\]
El parámetro \(\phi\) se especifica a través de la elección de la distribución, mientras que el parámetro de escala \(\sigma ^{2}\) permite una variabilidad extra. Por ejemplo, la Tabla 13.8 muestra que al especificar ya sea la distribución Bernoulli o Poisson, tenemos \(\phi =1\). Aunque el parámetro de escala \(\sigma ^{2}\) permite una variabilidad adicional, también puede acomodar situaciones en las que la variabilidad es menor de la especificada por la forma de la distribución (aunque esta situación es menos común). Finalmente, cabe destacar que para algunas distribuciones, como la normal, el término adicional ya está incorporado en el parámetro \(\phi\), por lo que no tiene un propósito útil.
Cuando se incluye el parámetro de escala adicional \(\sigma ^{2}\), es común estimarlo mediante el estadístico chi-cuadrado de Pearson dividido por los grados de libertad del error. Es decir,
\[
\widehat{\sigma }^{2}=\frac{1}{N-k} \sum_{i=1}^n w_i\frac{\left( y_i-b^{\prime }(\mathbf{x}_i^{\mathbf{\prime }}\mathbf{b}_{MLE})\right) ^{2}}{\phi b^{\prime \prime }(\mathbf{x}_i^{\mathbf{\prime }}\mathbf{b}_{MLE})}.
\]
Al igual que con la distribución Poisson en la Sección 12.3, otra forma de manejar patrones inusuales de varianza es a través de errores estándar robustos o empíricos. La Sección 13.9.3 proporciona más detalles.
Estadísticas de Bondad de Ajuste
En los modelos de regresión lineal, la estadística de bondad de ajuste más ampliamente citada es la medida \(R^2\), que se basa en la descomposición
\[
\sum_i \left(y_i - \overline{y} \right)^2 = \sum_i \left(y_i - \widehat{y}_i \right)^2 + \sum_i \left(\widehat{y}_i - \overline{y}\right)^2 + 2 \times \sum_i \left(y_i - \widehat{y}_i \right)\left(\widehat{y}_i - \overline{y}\right).
\]
En el lenguaje de la Sección 2.3, esta descomposición es:
\[
\textit{Total SS = Error SS + Regresión SS + 2} \times \textit{Suma de Productos Cruzados}.
\]
La dificultad con los modelos no lineales es que el término de Suma de Productos Cruzados rara vez es igual a cero. Por lo tanto, se obtienen estadísticas diferentes al definir \(R^2\) como (Regresión SS/Total SS) en comparación con (1-Error SS/Total SS). La Sección 11.3.2 describió algunas medidas alternativas de \(R^2\) que a veces se citan en configuraciones GLM.
Una medida de bondad de ajuste ampliamente citada es el estadístico chi-cuadrado de Pearson, que fue introducido en la Sección 12.1.4. En el contexto GLM, suponemos que \(\mathrm{E~}y_i = \mu_i\), \(\mathrm{Var~}y_i = \phi v(\mu_i)\) es la función de varianza (como en los ejemplos de la Tabla 13.1) y que \(\widehat{\mu}_i\) es un estimador de \(\mu_i\). Luego, el estadístico chi-cuadrado de Pearson se define como \(\sum_i \left( y_i - \widehat{\mu}_i \right)^2/\ ( \phi v(\widehat{\mu}_i))\). Como hemos visto en los modelos de Poisson para datos de conteo, esta formulación se reduce a la forma \(\sum_i \left( y_i - \widehat{\mu}_i \right)^2/\widehat{\mu}_i\).
Los criterios de información generales, incluyendo \(AIC\) y \(BIC\), que fueron definidos en la Sección 11.9 también se citan regularmente en estudios de GLM.
Una medida de bondad de ajuste específica para los modelos GLM es la estadística de devianza. Para definir esta estadística, trabajamos con la noción de un modelo saturado donde hay tantos parámetros como observaciones, \(\theta_i, i=1, \ldots, n\). Un modelo saturado proporciona el mejor ajuste posible. Con un parámetro para cada observación, maximizamos la verosimilitud en una base de observación por observación. Por lo tanto, al derivar la verosimilitud logarítmica de la ecuación (13.2) se obtiene
\[
\frac{\partial}{\partial \theta_i} \ln \mathrm{f}( y_i; \theta_i ,\phi ) = \frac{y_i-b^{\prime}(\theta_i)}{\phi}.
\]
Al igualar esto a cero, se obtiene la estimación del parámetro, digamos \(\theta_{i,SAT}\), como la solución de \(y_i=b^{\prime}(\theta_{i,Sat})\). Denotando \(\boldsymbol \theta _{SAT}\) como el vector de parámetros, la verosimilitud \(L(\boldsymbol \theta _{SAT})\) es el valor máximo posible de la verosimilitud logarítmica. Luego, para un estimador genérico \(\widehat{\boldsymbol \theta}\), la estadística de devianza escalada se define como
\[
\mathbf{\mathrm{D}}^{\ast}(\widehat{\boldsymbol \theta}) = 2 \times \left( L(\boldsymbol \theta _{SAT}) - L(\widehat{\boldsymbol \theta}) \right).
\]
En las familias exponenciales lineales, se multiplica por el factor de escala \(\phi\) para definir la estadística de devianza, \(\mathbf{\mathrm{D}}(\widehat{\boldsymbol \theta}) = \phi \mathbf{\mathrm{D}}^{\ast}(\widehat{\boldsymbol \theta})\). Esta multiplicación realmente elimina el factor de escala de la varianza de la definición de la estadística.
Es fácil verificar que la estadística de devianza se reduce a las siguientes formas para tres casos especiales:
- Normal: \(\mathbf{\mathrm{D}}(\widehat{\boldsymbol \mu})= \sum_i \left( y_i - \widehat \mu _i \right) ^2\),
- Bernoulli: \(\mathbf{\mathrm{D}}(\widehat{\boldsymbol \pi})= \sum_i \left\{ y_i \ln \frac {y_i}{\widehat \pi_i} + (1-y_i) \ln \frac {1-y_i}{1-\widehat \pi_i} \right\}\), y
- Poisson: \(\mathbf{\mathrm{D}}(\widehat{\boldsymbol \mu})= \sum_i \left\{ y_i \ln \frac {y_i}{\widehat \mu_i} + (y_i - \widehat \mu _i) \right\}\).
Aquí usamos la convención de que \(y \ln y = 0\) cuando \(y = 0\).
Aplicación: Gastos Médicos
Ahora volvemos a los datos de la Encuesta del Panel de Gastos Médicos (MEPS) introducidos en la Sección 11.4. En esa sección, buscamos desarrollar un modelo para entender el evento de una admisión hospitalaria como paciente interno. En esta sección, ahora deseamos modelar la cantidad del gasto. En terminología actuarial, la Sección 11.4 consideraba la “frecuencia”, mientras que esta sección involucra la “severidad”.
De las 2,000 observaciones seleccionadas aleatoriamente del año 2003 consideradas en la Sección 11.4, solo \(n=157\) fueron admitidos al hospital durante el año. La Tabla 13.5 resume los datos utilizando las mismas variables explicativas que en la Tabla 11.4. Por ejemplo, la Tabla 13.5 muestra que la muestra es 72% femenina, casi 76% blanca y más del 91% tiene seguro. La tabla también muestra relativamente pocos gastos por parte de asiáticos, nativos americanos y personas sin seguro en nuestra muestra.
La Tabla 13.5 también proporciona los gastos medianos por variable categórica. Esta tabla sugiere que el género, una mala autoevaluación de la salud física y los ingresos bajos o negativos pueden ser determinantes importantes de la cantidad de los gastos médicos.
Tabla 13.5. Gastos Medianos por Variable Explicativa Basado en una Muestra de \(n=157\) con Gastos Positivos
\[
\scriptsize{
\begin{array}{lllrr}
\hline
\text{Categoría} & \text{Variable} & \text{Descripción} & \text{Porcentaje} & \text{Mediana} \\
& & & \text{de datos} & \text{Gasto} \\ \hline
& COUNTIP & \text{Número de gastos (mediana: 1.0)} \\
\text{Demografía} & AGE & \text{Edad en años entre } \\
& & \ \ \ \text{ 18 a 65 (mediana: 41.0)} \\
& GENDER & \text{1 si es mujer} & 72.0 & 5,546 \\
& & \text{0 si es hombre} & 28.0 & 7,313 \\
\text{Etnicidad} & ASIAN & \text{1 si es asiático} & 2.6 & 4,003 \\
& BLACK & \text{1 si es afroamericano} & 19.8 & 6,100 \\
& NATIVE & \text{1 si es nativo} & 1.9 & 2,310 \\
& WHITE & \textit{Nivel de referencia} & 75.6 & 5,695 \\
\text{Región} & NORTHEAST & \text{1 si es del noreste} & 18.5 & 5,833 \\
& MIDWEST & \text{1 si es del medio oeste} & 21.7 & 7,999 \\
& SOUTH & \text{1 si es del sur} & 40.8 & 5,595 \\
& WEST & \textit{Nivel de referencia} & 19.1 & 4,297 \\
\hline \text{Educación} & COLLEGE
& \text{1 si tiene título universitario o superior} & 23.6 & 5,611 \\
& \text{HIGHSCHOOL} & \text{1 si tiene título de secundaria} & 43.3 & 5,907 \\
&& \text{Nivel de referencia es menor} & 33.1 & 5,338 \\
&& \ \ \ \text{ que título de secundaria} & \\ \hline
\text{Salud física} & POOR & \text{1 si es mala} & 17.2 & 10,447 \\
\ \ \text{autocalificada} & FAIR & \text{1 si es regular} & 10.2 & 5,228 \\
& GOOD & \text{1 si es buena } & 31.2 & 5,032 \\
& VGOOD & \text{1 si es muy buena} & 24.8 & 5,546 \\
&& \text{Nivel de referencia es excelente salud} & 16.6 & 5,277 \\
\text{Salud mental} & MPOOR & \text{1 si es mala o regular} & 15.9 & 6,583 \\
\ \ \text{autocalificada} & & \text{0 si es buena a excelente} & 84.1 & 5,599 \\
\text{Cualquier} & ANYLIMIT & \text{1 si tiene alguna limitación funcional/de actividad} &
41.4 & 7,826 \\
\ \ \text{limitación} & & \text{0 en caso contrario} & 58.6 & 4,746 \\
\hline \text{Ingreso} && \text{Nivel de referencia es ingreso alto} & 21.7 & 7,271 \\
\text{comparado con} & MINCOME & \text{1 si es ingreso medio} & 26.8 & 5,851 \\
\text{línea de pobreza} & LINCOME & \text{1 si es ingreso bajo} & 16.6 & 6,909 \\
& NPOOR & \text{1 si es casi pobre} & 7.0 & 5,546 \\
& POORNEG & \text{si es pobre/ingreso negativo} & 28.0 & 4,097 \\
\hline \text{Seguro} & INSURE & \text{1 si tiene cobertura de salud pública/privada}
& 91.1 & 5,943 \\
\ \ \text{médico} & & \ \ \text{en cualquier mes de 2003} & & \\
& & \text{0 si no tiene seguro médico en 2003} & 8.9 & 2,668 \\
\hline
Total & & & 100.0 & 5,695 \\ \hline
\end{array}
}
\]
knitr::kable(2, caption = "Silly. Create a table just to update the counter...")
Tabla 13.2: Silly. Create a table just to update the counter…
| 2 |
knitr::kable(2, caption = "Silly.")
knitr::kable(2, caption = "Silly. ")
Como se señaló en el prefacio, la mayor parte del código para el libro, escrito alrededor de 2008, estaba en SAS. El código ilustrativo proporcionado aquí, en R, presenta algunas ligeras discrepancias con el original.
Hexpend <- read.csv("CSVData/HealthExpend.csv", header=TRUE)
# Tabla 13.5
#Hexpend$POSEXP <- 1*(Hexpend$EXPENDIP>0)
HexpendPos <- subset(Hexpend, EXPENDIP > 0)
# CREAR UNA FUNCIÓN CORTA PARA AHORRAR TRABAJO
fun2 <- function(group_var, summary_vars){
# Convertir group_var y summary_vars a caracteres si no lo son
group_var <- as.character(substitute(group_var))
summary_vars <- as.character(substitute(summary_vars))
temp <- doBy::summaryBy(
formula = as.formula(paste(summary_vars, "~", group_var)),
data = HexpendPos,
FUN = function(x) { c(m = median(x), num = length(x)) } )
temp1 <- temp[,c(1,3,2)]
temp1[,2] <- round(100*temp1[,2]/length(HexpendPos$EXPENDIP), digits=1)
colnames(temp1) <- NULL
temp1[,1] <- as.character(temp1[,1])
return(as.matrix(temp1) )
}
var1 <- fun2(GENDER,EXPENDIP)
var2 <- fun2(RACE,EXPENDIP)
var3 <- fun2(REGION,EXPENDIP)
var4 <- fun2(EDUC,EXPENDIP)
var5 <- fun2(PHSTAT,EXPENDIP)
#var6 <- fun2(MPOOR,EXPENDIP)
var7 <- fun2(ANYLIMIT,EXPENDIP)
var8 <- fun2(INCOME,EXPENDIP)
var9 <- fun2(insure,EXPENDIP)
tableout <- rbind(var1, var2, var3, var4,
var5, var7, var8, var9)
TableGen1(TableData=tableout,
TextTitle='Gastos Medianos por Variable Explicativa Basado en una Muestra de $n=157$ con Gastos Positivos',
Align='crr', ColumnSpec=1:3, ColWidth = ColWidth4 )
Tabla 13.5: Gastos Medianos por Variable Explicativa Basado en una Muestra de \(n=157\) con Gastos Positivos
|
1
|
0
|
28
|
7312.695
|
|
2
|
1
|
72
|
5546.470
|
|
1
|
ASIAN
|
2.5
|
4002.735
|
|
2
|
BLACK
|
19.7
|
6100.200
|
|
3
|
NATIV
|
1.9
|
2310.320
|
|
4
|
OTHER
|
1.9
|
4050.650
|
|
5
|
WHITE
|
73.9
|
5725.325
|
|
1
|
MIDWEST
|
21.7
|
7998.660
|
|
2
|
NORTHEAST
|
18.5
|
5832.860
|
|
3
|
SOUTH
|
40.8
|
5595.295
|
|
4
|
WEST
|
19.1
|
4296.540
|
|
1
|
COLLEGE
|
23.6
|
5610.76
|
|
2
|
HIGHSCH
|
43.3
|
5907.33
|
|
3
|
LHIGHSC
|
33.1
|
5338.04
|
|
1
|
EXCE
|
16.6
|
5277.160
|
|
2
|
FAIR
|
10.2
|
5228.465
|
|
3
|
GOOD
|
31.2
|
5031.960
|
|
4
|
POOR
|
17.2
|
10447.130
|
|
5
|
VGOO
|
24.8
|
5546.470
|
|
1
|
0
|
58.6
|
4745.89
|
|
2
|
1
|
41.4
|
7825.50
|
|
1
|
HINCOME
|
21.7
|
7270.780
|
|
2
|
LINCOME
|
16.6
|
6908.720
|
|
3
|
MINCOME
|
26.8
|
5851.160
|
|
4
|
NPOOR
|
7.0
|
5546.470
|
|
5
|
POOR
|
28.0
|
4096.545
|
|
1
|
0
|
8.9
|
2668.37
|
|
2
|
1
|
91.1
|
5943.25
|
Tabla 13.5 utiliza medianas en lugar de medias debido a que la distribución de los gastos está sesgada hacia la derecha. Esto es evidente en la Figura 13.1, que proporciona un histograma suave (conocido como “estimación de densidad por kernel”, ver Sección 15.2) para los gastos hospitalarios. Para distribuciones sesgadas, la mediana a menudo brinda una mejor idea del centro de la distribución que la media. La distribución está aún más sesgada de lo que sugiere esta figura, ya que el gasto más grande (que es $607,800) se omite en la representación gráfica.
# ELIMINAR TEMPORALMENTE EL OBS # 58
Hexpend <- read.csv("CSVData/HealthExpend.csv", header=TRUE)
Hexpend58 <- subset(Hexpend, !(EXPENDIP>600000))
plot(density(Hexpend58$EXPENDIP), main="", xlab="Gastos")
#hist(Hexpend58$EXPENDIP, main="", xlab="Gastos")
Un modelo de regresión gamma con un enlace logarítmico se ajustó a los gastos hospitalarios utilizando todas las variables explicativas. Los resultados de este ajuste de modelo aparecen en la Tabla 13.6. Aquí vemos que muchos de los posibles determinantes importantes de los gastos médicos no son estadísticamente significativos. Esto es común en el análisis de gastos, donde las variables ayudan a predecir la frecuencia, aunque no son tan útiles para explicar la severidad.
Debido a la colinealidad, hemos visto en modelos lineales que tener demasiadas variables en un modelo ajustado puede llevar a la insignificancia estadística de variables importantes e incluso causar que los signos se inviertan. Para un modelo más simple, eliminamos las variables de Asian, Native American y uninsured, ya que representan un pequeño subconjunto de nuestra muestra. También utilizamos solo la variable POOR para el estado de salud autoinformado y solo POORNEG para los ingresos, esencialmente reduciendo estas variables categóricas a variables binarias. La Tabla 13.6, bajo el encabezado “Modelo Reducido”, muestra los resultados de este ajuste del modelo. Este modelo tiene casi el mismo estadístico de bondad de ajuste, AIC, lo que sugiere que es una alternativa razonable. Bajo el ajuste del modelo reducido, las variables COUNTIP (conteo de pacientes hospitalizados), AGE, COLLEGE y POORNEG son variables estadísticamente significativas. Para estas variables significativas, los signos son intuitivamente razonables. Por ejemplo, el coeficiente positivo asociado con COUNTIP significa que a medida que aumenta el número de visitas hospitalarias, los gastos totales aumentan, como era de esperar.
Para otra especificación alternativa del modelo, la Tabla 13.6 también muestra un ajuste de un modelo inverso gaussiano con un enlace logarítmico. A partir del estadístico AIC, vemos que este modelo no se ajusta tan bien como el modelo de regresión gamma. Todas las variables son estadísticamente insignificantes, lo que dificulta la interpretación de este modelo.
Tabla 13.6. Comparación de Modelos de Regresión Gamma e Inversa Gaussiana
\[
\scriptsize{
\begin{array}{l|rr|rr|rr}
\hline & \text{Gamma} & &\text{Gamma}&&\text{Inversa}&\text{Gaussiana} \\
\hline & \text{Modelo} &\text{Modelo}&
\text{Modelo} & \text{Modelo}&\text{Modelo}&\text{Modelo} \\
& \text{Completo} &\text{Reducido}&
\text{Completo} & \text{Reducido}&\text{Reducido}&\text{Reducido} \\ \hline
& \scriptsize{Parámetro} & & \scriptsize{Parámetro} & & \scriptsize{Parámetro} & \\
\text{Efecto} & \scriptsize{Estimación} & t\text{-valor} & \scriptsize{Estimación} & t\text{-valor} &
\scriptsize{Estimación} & t\text{-valor} \\ \hline
Intercepto & 6.891 & 13.080 & 7.859 & 17.951 & 6.544 & 3.024 \\
COUNTIP & 0.681 & 6.155 & 0.672 & 5.965 & 1.263 & 0.989 \\
AGE & 0.021 & 3.024 & 0.015 & 2.439 & 0.018 & 0.727 \\
GENDER & -0.228 & -1.263 & -0.118 & -0.648 & 0.363 & 0.482 \\
ASIAN & -0.506 & -1.029 & & & & \\
BLACK & -0.331 & -1.656 & -0.258 & -1.287 & -0.321 & -0.577 \\
NATIVE & -1.220 & -2.217 & & & & \\
NORTHEAST & -0.372 & -1.548 & -0.214 & -0.890 & 0.109 & 0.165 \\
MIDWEST & 0.255 & 1.062 & 0.448 & 1.888 & 0.399 & 0.654 \\
SOUTH & 0.010 & 0.047 & 0.108 & 0.516 & 0.164 & 0.319 \\
\hline
COLLEGE & -0.413 & -1.723 & -0.469 & -2.108 & -0.367 & -0.606 \\
HIGHSCHOOL & -0.155 & -0.827 & -0.210 & -1.138 & -0.039 & -0.078 \\
\hline
POOR & -0.003 & -0.010 & 0.167 & 0.706 & 0.167 & 0.258 \\
FAIR & -0.194 & -0.641 & & & & \\
GOOD & 0.041 & 0.183 & & & & \\
VGOOD & 0.000 & 0.000 & & & & \\
MNHPOOR & -0.396 & -1.634 & -0.314 & -1.337 & -0.378 & -0.642 \\
ANYLIMIT & 0.010 & 0.053 & 0.052 & 0.266 & 0.218 & 0.287 \\
\hline
MINCOME & 0.114 & 0.522 & & & & \\
LINCOME & 0.536 & 2.148 & & & & \\
NPOOR & 0.453 & 1.243 & & & & \\
POORNEG & -0.078 & -0.308 & -0.406 & -2.129 & -0.356 & -0.595 \\
INSURE & 0.794 & 3.068 & & & & \\
\hline
Escala & 1.409 & 9.779 & 1.280 & 9.854 & 0.026 & 17.720 \\
\hline
\scriptsize{Log-Verosimilitud} & -1,558.67 && -1,567.93 && -1,669.02 \\
AIC & 3,163.34 && 3,163.86 && 3,366.04 \\
\hline
\end{array}
}
\]
Como se señaló en el prefacio, la mayor parte del código para el libro, escrito alrededor de 2008, estaba en SAS. El código ilustrativo proporcionado aquí, en R, presenta algunas ligeras discrepancias con el original.
Hexpend <- read.csv("CSVData/HealthExpend.csv", header=TRUE)
# Tabla 13.6
#POSEXP = 1*(Hexpend$EXPENDIP>0)
HexpendPos <- subset(Hexpend, EXPENDIP > 0)
# MODELO CON TODAS LAS VARIABLES
# SE TUVO QUE AUMENTAR EL NÚMERO POR DEFECTO DE ITERACIONES PARA LA CONVERGENCIA
model1 = glm(EXPENDIP~COUNTIP+AGE+GENDER
+factor(RACE)+factor(REGION)+factor(EDUC)+factor(PHSTAT)
+MNHPOOR+ANYLIMIT+factor(INCOME)+insure,
control = glm.control(maxit = 50),
data=HexpendPos,family=Gamma(link="log"))
output1 <- summary(model1 )
tidy_model <- broom::tidy(model1)
kable(tidy_model, format = "html", digits = 3,
table.attr = "class='table table-striped'",
caption = "Resumen del Modelo de Regresión Gamma")
Tabla 13.6: Resumen del Modelo de Regresión Gamma
|
term
|
estimate
|
std.error
|
statistic
|
p.value
|
|
(Intercept)
|
6.211
|
0.664
|
9.352
|
0.000
|
|
COUNTIP
|
0.682
|
0.094
|
7.287
|
0.000
|
|
AGE
|
0.021
|
0.007
|
2.989
|
0.003
|
|
GENDER
|
-0.238
|
0.172
|
-1.382
|
0.169
|
|
factor(RACE)BLACK
|
0.182
|
0.512
|
0.356
|
0.722
|
|
factor(RACE)NATIV
|
-0.697
|
0.732
|
-0.953
|
0.343
|
|
factor(RACE)OTHER
|
0.726
|
0.707
|
1.027
|
0.306
|
|
factor(RACE)WHITE
|
0.513
|
0.481
|
1.066
|
0.289
|
|
factor(REGION)NORTHEAST
|
-0.628
|
0.238
|
-2.642
|
0.009
|
|
factor(REGION)SOUTH
|
-0.256
|
0.199
|
-1.288
|
0.200
|
|
factor(REGION)WEST
|
-0.261
|
0.233
|
-1.118
|
0.266
|
|
factor(EDUC)HIGHSCH
|
0.261
|
0.193
|
1.356
|
0.177
|
|
factor(EDUC)LHIGHSC
|
0.411
|
0.225
|
1.828
|
0.070
|
|
factor(PHSTAT)FAIR
|
-0.165
|
0.314
|
-0.524
|
0.601
|
|
factor(PHSTAT)GOOD
|
0.063
|
0.227
|
0.276
|
0.783
|
|
factor(PHSTAT)POOR
|
0.026
|
0.298
|
0.088
|
0.930
|
|
factor(PHSTAT)VGOO
|
0.018
|
0.237
|
0.076
|
0.939
|
|
MNHPOOR
|
-0.397
|
0.238
|
-1.664
|
0.098
|
|
ANYLIMIT
|
-0.006
|
0.189
|
-0.031
|
0.976
|
|
factor(INCOME)LINCOME
|
0.551
|
0.257
|
2.147
|
0.034
|
|
factor(INCOME)MINCOME
|
0.123
|
0.211
|
0.583
|
0.561
|
|
factor(INCOME)NPOOR
|
0.463
|
0.354
|
1.308
|
0.193
|
|
factor(INCOME)POOR
|
-0.068
|
0.251
|
-0.273
|
0.785
|
|
insure
|
0.789
|
0.257
|
3.068
|
0.003
|
# MODELO REDUCIDO
HexpendPos$BLACK = (HexpendPos$RACE == "BLACK")
HexpendPos$POOR = (HexpendPos$PHSTAT == "POOR")
HexpendPos$POORNEG = (HexpendPos$INCOME == "POOR")
model2 <- glm(EXPENDIP ~ COUNTIP + AGE + GENDER +
factor(BLACK) + factor(REGION) + factor(EDUC) + factor(POOR) +
MNHPOOR + ANYLIMIT + factor(POORNEG),
data=HexpendPos, family=Gamma(link="log") )
#summary(model2)
anova(model1, model2)
Analysis of Deviance Table
Model 1: EXPENDIP ~ COUNTIP + AGE + GENDER + factor(RACE) + factor(REGION) +
factor(EDUC) + factor(PHSTAT) + MNHPOOR + ANYLIMIT + factor(INCOME) +
insure
Model 2: EXPENDIP ~ COUNTIP + AGE + GENDER + factor(BLACK) + factor(REGION) +
factor(EDUC) + factor(POOR) + MNHPOOR + ANYLIMIT + factor(POORNEG)
Resid. Df Resid. Dev Df Deviance
1 133 123.954
2 143 137.876 -10 -13.9218
Hexpend58 <- subset(HexpendPos, !(EXPENDIP>600000))
# MODELO REDUCIDO - INVERSA GAUSSIANA
# ¡NO QUIERE CONVERGER!
# model3 <- glm(EXPENDIP ~ COUNTIP + AGE + GENDER +
# factor(BLACK) + factor(REGION) + factor(EDUC) + factor(POOR) +
# MNHPOOR + ANYLIMIT + factor(POORNEG),
# data=Hexpend58, family=inverse.gaussian(link = "log") )
# summary(model3)
# TAMBIÉN PROBLEMAS CON family=inverse.gaussian(link = "1/mu^2")
Residuales
Una forma de seleccionar un modelo adecuado es ajustar una especificación tentativa y analizar formas de mejorarlo. Este método diagnóstico se basa en los residuales. En el modelo lineal, definimos los residuales como la respuesta menos los valores ajustados correspondientes. En un contexto de GLM, nos referimos a estos como “residuales crudos”, denotados como \(y_i - \widehat{\mu}_i\). Como hemos visto, los residuales son útiles para:
- descubrir nuevas covariables o los efectos de patrones no lineales en covariables existentes,
- identificar observaciones con mal ajuste,
- cuantificar los efectos de las observaciones individuales en los parámetros del modelo, y
- revelar otros patrones del modelo, como heterocedasticidad o tendencias temporales.
Como se enfatiza en esta sección, los residuales también son los componentes básicos de muchas estadísticas de bondad de ajuste.
Como vimos en la Sección 11.1 sobre variables dependientes binarias, el análisis de los residuales crudos puede ser irrelevante en algunos contextos no lineales. Cox y Snell (1968, 1971) introdujeron una noción general de residuales que es más útil en modelos no lineales que los simples residuales crudos. Para nuestras aplicaciones, asumimos que \(y_i\) tiene una función de distribución F(\(\cdot\)) que está indexada por las variables explicativas \(\mathbf{x}_i\) y un vector de parámetros \(\boldsymbol \theta\) denotado como F(\(\mathbf{x}_i,\boldsymbol \theta\)). Aquí, la función de distribución es común a todas las observaciones, pero la distribución varía a través de las variables explicativas. El vector de parámetros \(\boldsymbol \theta\) incluye los parámetros de regresión \(\boldsymbol \beta\) así como los parámetros de escala. Con conocimiento de F(\(\cdot\)), ahora asumimos que se puede calcular una nueva función R(\(\cdot\)) tal que \(\varepsilon_i = \mathrm{R}(y_i; \mathbf{x}_i,\boldsymbol \theta)\), donde \(\varepsilon_i\) son distribuidos idénticamente e independientemente. Aquí, la función R(\(\cdot\)) depende de las variables explicativas \(\mathbf{x}_i\) y los parámetros \(\boldsymbol \theta\). Los residuales de Cox-Snell son entonces \(e_i = \mathrm{R}(y_i; \mathbf{x}_i,\widehat{\boldsymbol \theta})\). Si el modelo está bien especificado y las estimaciones de los parámetros \(\widehat{\boldsymbol \theta}\) están cerca de los parámetros del modelo \(\boldsymbol \theta\), entonces los residuales \(e_i\) deberían estar cerca de ser i.i.d.
Para ilustrar, en el modelo lineal utilizamos la función de residuales
\(\mathrm{R}(y_i; \mathbf{x}_i,\boldsymbol \theta)\) \(= y_i - \mathbf{x}_i^{\prime} \boldsymbol \beta\) \(= y_i - \mu_i = \varepsilon_i\), resultando en residuales crudos (ordinarios). Otra opción es reescalar mediante la desviación estándar
\[
\mathrm{R}(y_i; \mathbf{x}_i,\boldsymbol \theta) = \frac{y_i -
\mu_i}{\sqrt{\mathrm{Var~} y_i}} ,
\]
que produce residuales de Pearson.
Otro tipo de residual se basa en transformar primero la respuesta con una función conocida h(\(\cdot\)),
\[
\mathrm{R}(y_i; \mathbf{x}_i,\boldsymbol \theta) =
\frac{\mathrm{h}(y_i) - \mathrm{E~}
\mathrm{h}(y_i)}{\sqrt{\mathrm{Var~}\mathrm{h}(y_i)}}.
\]
Elegir la función de modo que h(\(y_i\)) sea aproximadamente distribuida normalmente produce residuales de Anscombe. Otra elección popular es seleccionar la transformación de modo que la varianza sea estabilizada. Tabla 13.7 presenta transformaciones para las distribuciones binomial, Poisson y gamma.
Tabla 13.7. Transformaciones Aproximadas a la Normalidad y de Estabilización de Varianza
\[
\small{
\begin{array}{lcc}
\hline
\text{Distribución} & \text{Transformación Aproximada } & \text{Transformación de Estabilización } \\
& \text{a la Normalidad Transformación} & \text{de Varianza Transformación} \\
\hline \text{Binomial}(m,p) &
\frac{\mathrm{h}_B(y/m)-[\mathrm{h}_B(p)+(p(1-p))^{-1/3}(2p-1)/(6m)]}{(p(1-p)^{1/6}/\sqrt{m})}
&
\frac{\sin^{-1}(\sqrt{y/m})-\sin^{-1}(\sqrt{p})}{a/(2\sqrt{m})} \\
\text{Poisson}(\lambda) & 1.5 \lambda^{-1/6} \left[
y^{2/3}-(\lambda^{2/3}-\lambda^{-1/3}/9) \right] &
2(\sqrt{y}-\sqrt{\lambda})\\
\text{Gamma}(\alpha, \gamma) & 3 \alpha^{1/6} \left[
(y/\gamma)^{1/3}-(\alpha^{1/3}-\alpha^{-2/3}/9) \right]
& \text{No considerado}\\
\hline \\
\end{array}
}
\]
Nota: \(\mathrm{h}_B(u) = \int_0^u s^{-1/3} (1-s)^{-1/3}ds\) es la función beta incompleta evaluada en \((2/3, 2/3)\), hasta una constante escalar. Fuente: Pierce y Schafer (1986).
Otra opción ampliamente utilizada en estudios de GLM son los residuales de devianza, definidos como
\[
\mathrm{sign} \left(y_i -
\widehat{\mu}_i \right)
\sqrt{2 \left( \ln \mathrm{f}(y_i;\theta_{i,SAT}) - \ln
\mathrm{f}(y_i;\widehat{\theta}_{i})
\right)}.
\]
La Sección 13.3.3 introdujo las estimaciones de los parámetros del modelo saturado, \(\theta_{i,SAT}\). Como mencionan Pierce y Schafer (1986), los residuales de devianza son muy cercanos a los residuales de Anscombe en muchos casos; podemos verificar los residuales de devianza para la normalidad aproximada. Además, los residuales de devianza pueden definirse fácilmente para cualquier modelo GLM y son fáciles de calcular.
Distribución de Tweedie
Hemos visto que la familia exponencial natural incluye distribuciones continuas, como la normal y gamma, así como distribuciones discretas, como la binomial y Poisson. También incluye distribuciones que son mezclas de componentes discretos y continuos. En la modelización de reclamaciones de seguros, la mezcla más utilizada es la distribución de Tweedie (1984). Esta tiene una masa positiva en cero, representando la ausencia de reclamaciones, y un componente continuo para valores positivos que representan el monto de una reclamación.
La distribución de Tweedie se define como una suma de Poisson de variables aleatorias gamma, conocida como una pérdida agregada en ciencias actuariales. Específicamente, supongamos que \(N\) tiene una distribución de Poisson con media \(\lambda\), representando el número de reclamaciones. Sea {\(y_j\)} una secuencia i.i.d., independiente de \(N\), con cada \(y_j\) siguiendo una distribución gamma con parámetros \(\alpha\) y \(\gamma\), representando el monto de una reclamación. Entonces, \(S_N = y_1 + \ldots + y_N\) es una suma de Poisson de gammas. La Sección 16.5 discutirá con más detalle los modelos de pérdidas agregadas. Esta sección se enfoca en el caso especial importante de la distribución de Tweedie.
Para entender el aspecto de mezcla de la distribución de Tweedie, primero observa que es sencillo calcular la probabilidad de tener cero reclamaciones como
\[
\Pr ( S_N=0) = \Pr (N=0) = e^{-\lambda}.
\]
La función de distribución puede calcularse usando expectativas condicionales,
\[
\Pr ( S_N \le y) = e^{-\lambda}+\sum_{n=1}^{\infty} \Pr(N=n) \Pr(S_n \le y), ~~~~~y \geq 0.
\]
Dado que la suma de gammas i.i.d. es una gamma, \(S_n\) (no \(S_N\)) tiene una distribución gamma con parámetros \(n\alpha\) y \(\gamma\). Así, para \(y>0\), la densidad de la distribución de Tweedie es
\[\begin{equation}
\mathrm{f}_{S}(y) = \sum_{n=1}^{\infty} e^{-\lambda} \frac{\lambda ^n}{n!} \frac{\gamma^{n \alpha}}{\Gamma(n \alpha)} y^{n \alpha -1} e^{-y \gamma}.
\tag{13.9}
\end{equation}\]
A primera vista, esta densidad no parece pertenecer a la familia exponencial lineal dada en la ecuación (13.2). Para ver la relación, primero calculamos los momentos usando expectativas iteradas como
\[\begin{equation}
\mathrm{E~}S_N = \lambda \frac{\alpha}{\gamma}~~~~~~\mathrm{y}~~~~~ \mathrm{Var~}S_N = \frac{\lambda \alpha}{\gamma^2} (1+\alpha).
\tag{13.10}
\end{equation}\]
Ahora, definimos tres parámetros \(\mu, \phi, p\) mediante las relaciones
\[
\lambda = \frac{\mu^{2-p}}{\phi (2-p)},~~~~~~\alpha = \frac{2-p}{p-1}~~~~~~\mathrm{y}~~~~~ \frac{1}{\gamma} = \phi(p-1)\mu^{p-1}.
\]
Al insertar estos nuevos parámetros en la ecuación (13.9) obtenemos
\[\begin{equation}
\mathrm{f}_{S}(y) = \exp \left[ \frac{-1}{\phi} \left( \frac{\mu^{2-p}}{2-p} + \frac{y}{(p-1)\mu^{p-1}} \right) + S(y,\phi) \right].
\tag{13.11}
\end{equation}\]
Dejamos el cálculo de \(S(y,\phi)\) como un ejercicio.
Así, la distribución de Tweedie pertenece a la familia exponencial lineal. Cálculos sencillos muestran que
\[\begin{equation}
\mathrm{E~}S_N = \mu~~~~~~\mathrm{y}~~~~~ \mathrm{Var~}S_N = \phi \mu^p,
\tag{13.12}
\end{equation}\]
donde \(1<p<2.\) Al examinar nuestra función de varianza Tabla 13.1, la distribución de Tweedie también puede verse como una elección que es intermedia entre las distribuciones Poisson y gamma.
Lecturas adicionales y referencias
Un tratamiento más extenso de los GLM puede encontrarse en el trabajo clásico de McCullagh y Nelder (1989) o la introducción más sencilla de Dobson (2002). de Jong y Heller (2008) proporcionan una reciente introducción a los GLM con un énfasis especial en seguros.
Los GLM han recibido una considerable atención de parte de los actuarios de seguros de propiedad y accidentes (generales) en los últimos tiempos. Véase Anderson et al. (2004) y Clark y Thayer (2004) para introducciones.
Mildenhall (1999) proporciona un análisis sistemático de la conexión entre los modelos GLM y los algoritmos para evaluar diferentes planes de tarificación. Fu y Wu (2008) proporcionan una versión actualizada, introduciendo una clase más amplia de algoritmos iterativos que se conectan bien con los modelos GLM.
Para más información sobre la distribución de Tweedie, consulte J{}rgensen y de Souza (1994) y Smyth y J{}rgensen (2002).
Referencias del capítulo
- Anderson, Duncan, Sholom Feldblum, Claudine Modlin, Doris
Schirmacher, Ernesto Schirmacher and Neeza Thandi (2004). A
practitioner’s guide to generalized linear models. Casualty Actuarial Society 2004 Discussion Papers 1-116, Casualty Actuarial Society.
- Baxter, L. A., S. M. Coutts and G. A. F. Ross (1980). Applications of linear models in motor insurance. Proceedings of the 21st International Congress of Actuaries. Zurich. pp. 11-29.
- Clark, David R. and Charles A. Thayer (2004). A primer on the exponential family of distributions. Casualty Actuarial Society 2004 Discussion Papers 117-148, Casualty Actuarial Society.
- Cox, David R. and E. J. Snell (1968). A general definition of
residuals. Journal of the Royal Statistical Society, Series
B, 30, 248-275.
- Cox, David R. and E. J. Snell (1971). On test statistics computed
from residuals. Biometrika 71, 589-594.
- Dobson, Annette J. (2002). An Introduction to Generalized Linear Models, Second Edition. Chapman & Hall, London.
- Fu, Luyang and Cheng-sheng Peter Wu (2008). General iteration algorithm for classification ratemaking. Variance 1(2), 193-213, Casualty Actuarial Society.
- de Jong, Piet and Gillian Z. Heller (2008). Generalized Linear Models for Insurance Data. Cambridge University Press, Cambridge, UK.
- J{}rgensen, Bent and Marta C. Paes de Souza (1994). Fitting Tweedie’s compound Poisson model to insurance claims data. Scandinavian Actuarial Journal 1994;1, 69-93.
- McCullagh, P. and J.A. Nelder (1989). Generalized Linear Models, Second Edition. Chapman & Hall, London.
- Mildenhall, Stephen J. (1999). A systematic relationship between minimum bias and generalized linear models. Casualty Actuarial Society Proceedings 86, Casualty Actuarial Society.
- Pierce, Donald A. and Daniel W. Schafer (1986). Residuals in generalized linear models. Journal of the American Statistical Association 81, 977-986.
- Smyth, Gordon K. and Bent J{}rgensen (2002). Fitting Tweedie’s compound Poisson model to insurance claims data: Dispersion modelling. Astin Bulletin 32(1), 143-157.
- Tweedie, M. C. K. (1984). An index which distinguishes between some important exponential families. In Statistics: Applications and New Directions. Proceedings of the Indian Statistical Golden Jubilee International Conference* (Editors J. K. Ghosh and J. Roy), pp. 579-604. Indian Statistical Institute, Calcutta.
Ejercicios
13.1 Verifica las entradas en Tabla 13.8 para la distribución gamma. Específicamente:
Demuestra que la distribución gamma es un miembro de la familia exponencial lineal de distribuciones.
Describe los componentes de la familia exponencial lineal (\(\theta, \phi, b(\theta), S(y,\phi)\)) en términos de los parámetros de la distribución gamma.
Demuestra que la media y la varianza de la distribución gamma y de las distribuciones de la familia exponencial lineal coinciden.
13.2 Clasificación de Tarificación. Este ejercicio considera los datos descritos en el ejemplo de clasificación de tarificación de la Sección 13.2.2, utilizando los datos de la Tabla 13.3.
Ajusta un modelo de regresión gamma usando una función de enlace logarítmica con los conteos de reclamaciones como ponderaciones (\(w_i\)). Utiliza las variables explicativas categóricas de grupo de edad y uso del vehículo para estimar la severidad media esperada.
Con base en tus parámetros estimados en la parte (a), verifica las severidades esperadas estimadas en la Tabla 13.4.
13.3 Verifica que la distribución de Tweedie es un miembro de la familia exponencial lineal de distribuciones comprobando la ecuación (13.9). En particular, proporciona una expresión para \(S(y,\phi)\) (nota que \(S(y,\phi)\) también depende de \(p\) pero no de \(\mu\)). Puedes consultar a Clark y Thayer (2004) para verificar tu trabajo. Además, verifica los momentos en la ecuación (13.10).
Suplementos Técnicos - Familia Exponencial
Familia Exponencial Lineal de Distribuciones
La distribución de la variable aleatoria \(y\) puede ser discreta, continua o una mezcla. Así, f(.) en la ecuación (13.2) puede interpretarse como una función de densidad o una función de masa, dependiendo de la aplicación. Tabla 13.8 proporciona varios ejemplos, incluyendo las distribuciones normal, binomial y Poisson.
Tabla 13.8. Distribuciones Seleccionadas de la Familia Exponencial de Un Parámetro
\[
\scriptsize{
\begin{array}{l|ccccc}
\hline
& & \text{Densidad o} & & & \\
\text{Distribución} & \text{Parámetros} & \text{Función de Masa} & \text{Componentes} & \mathrm{E}~y & \mathrm{Var}~y \\
\hline \text{General} & \theta,~ \phi &
\exp \left( \frac{y\theta -b(\theta )}{\phi } +S\left( y,\phi \right) \right)
& \theta,~ \phi, b(\theta), S(y, \phi) & b^{\prime}(\theta) & b^{\prime \prime}(\theta) \phi \\
\text{Normal} & \mu, \sigma^2 &
\frac{1}{\sigma \sqrt{2\pi }}\exp \left( -\frac{(y-\mu )^{2}}{2\sigma ^{2}}\right) & \mu, \sigma^2,
\frac{\theta^2}{2}, - \left(\frac{y^2}{2\phi} + \frac{\ln(2 \pi \phi)}{2} \right) & \theta=\mu & \phi= \sigma^2 \\
\text{Binomial} & \pi & {n \choose y} \pi ^y (1-\pi)^{n-y} & \ln
\left(\frac{\pi}{1-\pi} \right), 1, n \ln(1+e^{\theta} ), & n
\frac {e^{\theta}}{1+e^{\theta}} & n \frac
{e^{\theta}}{(1+e^{\theta})^2} \\
& & & \ln {n \choose y} & = n \pi & = n \pi (1-\pi) \\
\text{Poisson} & \lambda & \frac{\lambda^y}{y!} \exp(-\lambda) &
\ln \lambda, 1, e^{\theta}, - \ln (y!) & e^{\theta} = \lambda &
e^{\theta} = \lambda \\
\text{Negativa} & r,p & \frac{\Gamma(y+r)}{y! \Gamma(r)} p^r ( 1-p)^y & \ln(1-p), 1, -r \ln(1-e^{\theta}), &
\frac{r(1-p)}{p}
& \frac{r(1-p)}{p^2} \\
\ \ \ \text{Binomial}^{\ast} & & & ~~~\ln \left[ \frac{\Gamma(y+r)}{y! \Gamma(r)} \right] & = \mu &
= \mu+\mu^2/r
\\
\text{Gamma} & \alpha, \gamma & \frac{\gamma ^ \alpha}{\Gamma (\alpha)}
y^{\alpha -1 }\exp(-y \gamma) & - \frac{\gamma}{\alpha},
\frac{1}{\alpha}, - \ln ( - \theta), -\phi^{-1} \ln \phi & -
\frac{1}{\theta} = \frac{\alpha}{\gamma} & \frac{\phi}{\theta ^2}
=
\frac{\alpha}{\gamma ^2} \\
& & & - \ln \left( \Gamma(\phi ^{-1}) \right) +
(\phi^{-1} - 1) \ln y & & \\
\text{Inversa} & \mu, \lambda & \sqrt{\frac{\lambda}{2 \pi y^3}
} \exp \left(- \frac{\lambda (y-\mu)^2}{2 \mu^2 y} \right) &
-1/( 2\mu^2), 1/\lambda, -\sqrt{-2 \theta}, & \left(-2 \theta
\right)^{-1/2} & \phi (-2 \theta)^{-3/2} \\
\ \ \ \text{Gaussiana} & & & \theta/(\phi y) - 0.5 \ln (\phi 2 \pi y^3 ) & = \mu & = \frac{\mu^3}{\lambda} \\
\text{Tweedie} & & \textit{Ver Sección 13.6}\\
\hline \\
\end{array}
}
\]
\(^{\ast}\) Esto supone que el parámetro \(r\) es fijo, pero no tiene que ser un número entero.
Momentos
Para evaluar los momentos de las familias exponenciales, es conveniente trabajar con la función generadora de momentos. Para simplificar, asumimos que la variable aleatoria \(y\) es continua. Definimos la función generadora de momentos
\[\begin{eqnarray*}
M(s) &=& \mathrm{E}~ e^{sy} = \int \exp \left( sy+ \frac{y \theta - b(\theta)}{\phi} + S(y, \phi) \right) dy \\
&=& \exp \left(\frac{b(\theta + s \phi) - b(\theta)}{\phi} \right) \int \exp \left( \frac{y(\theta + s\phi) - b(\theta+s\phi)}{\phi} + S(y, \phi) \right) dy \\
&=& \exp \left(\frac{b(\theta + s \phi) - b(\theta)}{\phi} \right),
\end{eqnarray*}\]
porque \(\int \exp \left( \frac{1}{\phi} \left[ y(\theta + s\phi) - b(\theta+s\phi) \right] + S(y, \phi) \right) dy=1\). Con esta expresión, podemos generar los momentos. Así, para la media, tenemos
\[\begin{eqnarray*}
\mathrm{E}~ y &=& M^{\prime}(0) = \left . \frac{\partial}{\partial s} \exp \left(\frac{b(\theta + s \phi) - b(\theta)}{\phi} \right) \right | _{s=0} \\
&=& \left[ b^{\prime}(\theta + s \phi ) \exp \left( \frac{b(\theta + s \phi) - b(\theta)}{\phi} \right) \right] _{s=0} = b^{\prime}(\theta).
\end{eqnarray*}\]
De manera similar, para el segundo momento, tenemos
\[\begin{eqnarray*}
M^{\prime \prime}(s) &=& \frac{\partial}{\partial s} \left[ b^{\prime}(\theta + s \phi) \exp \left(\frac{b(\theta + s \phi) - b(\theta)}{\phi} \right) \right] \\
&=& \phi b^{\prime \prime}(\theta + s \phi) \exp \left(\frac{b(\theta + s \phi) - b(\theta)}{\phi} \right) + (b^{\prime}(\theta + s \phi))^2 \left(\frac{b(\theta + s \phi) - b(\theta)}{\phi} \right).
\end{eqnarray*}\]
Esto produce \(\mathrm{E}y^2 = M^{\prime \prime}(0)\) \(= \phi b^{\prime \prime}(\theta) + (b^{\prime}(\theta))^2\) y \(\mathrm{Var}~y =\phi b^{\prime \prime}(\theta)\).
Estimación de Máxima Verosimilitud para Enlaces Generales
Para enlaces generales, ya no asumimos la relación \(\theta_i=\mathbf{x}_i^{\mathbf{\prime }}\boldsymbol \beta\), sino que asumimos que \(\boldsymbol \beta\) está relacionado con \(\theta_i\) a través de las relaciones \(\mu _i=b^{\prime }(\theta_i)\) y \(\mathbf{x}_i^{\mathbf{\prime }}\boldsymbol \beta = g\left( \mu _i\right)\). Seguimos asumiendo que el parámetro de escala varía según la observación, de modo que \(\phi_i = \phi/w_i\), donde \(w_i\) es una función de ponderación conocida. Usando la ecuación (13.4), tenemos que el \(j\)-ésimo elemento de la función de puntaje es
\[
\frac{\partial }{\partial \beta _{j}}\ln \mathrm{f}\left( \mathbf{y}\right) =\sum_{i=1}^n \frac{\partial \theta_i}{\partial \beta _{j}} \frac{y_i-\mu _i}{\phi _i},
\]
ya que \(b^{\prime }(\theta_i)=\mu _i\). Ahora, usamos la regla de la cadena y la relación \(\mathrm{Var~}y_i=\phi _ib^{\prime \prime }(\theta_i)\) para obtener
\[
\frac{\partial \mu _i}{\partial \beta _{j}}=\frac{\partial b^{\prime }(\theta_i)}{\partial \beta _{j}}=b^{\prime \prime }(\theta_i)\frac{ \partial \theta_i}{\partial \beta _{j}}=\frac{\mathrm{Var~}y_i}{\phi _i}\frac{\partial \theta_i}{\partial \beta _{j}}.
\]
Así, tenemos
\[
\frac{\partial \theta_i}{\partial \beta _{j}}\frac{1}{\phi _i}=\frac{\partial \mu _i}{\partial \beta _{j}}\frac{1}{\mathrm{Var~}y_i}.
\]
Esto produce
\[
\frac{\partial }{\partial \beta _{j}}\ln \mathrm{f}\left( \mathbf{y}\right) =\sum_{i=1}^n \frac{\partial \mu _i}{\partial \beta _{j}}\left( \mathrm{Var~}y_i\right) ^{-1}\left( y_i-\mu _i\right),
\]
que resumimos como
\[\begin{equation}
\frac{\partial L(\boldsymbol \beta) }{\partial \boldsymbol \beta } = \frac{\partial }{\partial \boldsymbol \beta }\ln \mathrm{f}\left( \mathbf{y}\right) =\sum_{i=1}^n \frac{\partial \mu _i}{\partial \boldsymbol \beta}\left( \mathrm{Var~}y_i\right) ^{-1}\left( y_i-\mu _i\right),
\tag{13.13}
\end{equation}\]
que se conoce como la forma de las ecuaciones de estimación generalizadas.
Resolver las raíces de la ecuación de puntaje \(L(\boldsymbol \beta) = \mathbf{0}\) proporciona estimaciones de máxima verosimilitud, \(\mathbf{b}_{MLE}\). En general, esto requiere métodos numéricos iterativos. Una excepción es el siguiente caso especial.
Caso especial. Sin covariables. Supongamos que \(k=1\) y que \(x_i =1\). Entonces, \(\beta = \eta = g(\mu)\), y \(\mu\) no depende de \(i\). A partir de la relación \(\mathrm{Var~}y_i = \phi b^{\prime \prime}(\theta)/w_i\) y la ecuación (13.13), la función de puntaje es
\[
\frac{\partial L( \beta) }{\partial \beta } = \frac{\partial \mu}{\partial \beta} ~ \frac{1}{\phi b^{\prime \prime}(\theta)} \sum_{i=1}^n w_i \left( y_i-\mu \right).
\]
Igualando esto a cero se obtiene \(\overline{y}_w = \sum_i w_i y_i / \sum_i w_i = \widehat{\mu}_{MLE} = g^{-1}(b_{MLE})\), o \(b_{MLE}=g(\overline{y}_w)\), donde \(\overline{y}_w\) es un promedio ponderado de los \(y\).
Para la matriz de información, usamos la independencia entre los sujetos y la ecuación (13.13) para obtener
\[\begin{eqnarray}
\mathbf{I}(\boldsymbol \beta)
&=& \mathrm{E~} \left( \frac{\partial L(\boldsymbol \beta) }{\partial \boldsymbol \beta } \frac{\partial L(\boldsymbol \beta) }{\partial \boldsymbol \beta ^{\prime}} \right) \nonumber \\
&=& \sum_{i=1}^n \frac{\partial \mu _i}{\partial \boldsymbol \beta}~\mathrm{E} \left( \left( \mathrm{Var~}y_i\right) ^{-1}\left( y_i-\mu _i\right)^2 \left( \mathrm{Var~}y_i\right) ^{-1}\right) \frac{\partial \mu _i}{\partial \boldsymbol \beta ^{\prime}} \nonumber \\
&=& \sum_{i=1}^n \frac{\partial \mu _i}{\partial \boldsymbol \beta} \left( \mathrm{Var~}y_i\right) ^{-1} \frac{\partial \mu _i}{\partial \boldsymbol \beta ^{\prime}} .
\tag{13.14}
\end{eqnarray}\]
Tomar la raíz cuadrada de los elementos diagonales de la inversa de \(\mathbf{I}(\boldsymbol \beta)\), después de insertar las estimaciones de los parámetros en las funciones de media y varianza, proporciona errores estándar basados en el modelo. Para datos en los que se duda de la correcta especificación del componente de varianza, se reemplaza Var \(y_i\) por una estimación insesgada \((y_i - \mu_i)^2\). El estimador resultante
\[
se(b_j) = \sqrt{ j\text{-ésimo elemento diagonal de } [\sum_{i=1}^n \frac{\partial \mu _i}{\partial \boldsymbol \beta} (y_i - \mu_i)^{-2} \frac{\partial \mu _i}{\partial \boldsymbol \beta ^{\prime}}]|_{\boldsymbol \beta = \mathbf{b}_{MLE}}^{-1}}
\]
se conoce como el error estándar empírico o robusto.
Mínimos Cuadrados Reponderados Iterativos
Para ver cómo funcionan los métodos iterativos generales, usamos el enlace canónico de modo que \(\eta_i = \theta_i\). Luego, la función de puntaje está dada en la ecuación (13.6) y el Hessiano en la ecuación (13.8). Usando la ecuación de Newton-Raphson (11.15) en la Sección 11.9 obtenemos
\[\begin{equation}
\boldsymbol \beta_{NEW} = \boldsymbol \beta_{OLD} - \left( \sum_{i=1}^n w_i b^{\prime \prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta _{OLD}) \mathbf{x}_i \mathbf{x}_i^{\prime} \right)^{-1} \left( \sum_{i=1}^n w_i (y_i - b^{\prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta_{OLD})) \mathbf{x}_i \right).
\tag{13.15}
\end{equation}\]
Observa que, dado que la matriz de segundas derivadas no es estocástica, Newton-Raphson es equivalente al algoritmo de puntuación de Fisher. Como se describe en la Sección 13.3.1, la estimación de \(\boldsymbol \beta\) no requiere el conocimiento de \(\phi\).
Seguimos usando el enlace canónico y definimos una “variable dependiente ajustada”
\[
y_i^{\ast}(\boldsymbol \beta) = \mathbf{x}_i^{\prime} \boldsymbol \beta + \frac{y_i - b^{\prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta)}{b^{\prime \prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta)}.
\]
Esto tiene varianza
\[
\mathrm{Var~}y_i^{\ast}(\boldsymbol \beta) = \frac{\mathrm{Var~}y_i}{(b^{\prime \prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta))^2} = \frac{\phi_i b^{\prime \prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta)}{(b^{\prime \prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta))^2} = \frac{\phi/ w_i }{b^{\prime \prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta)}.
\]
Usamos el nuevo peso como el recíproco de la varianza, \(w_i(\boldsymbol \beta)=w_i b^{\prime \prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta) / \phi.\) Luego, con la expresión
\[
w_i \left(y_i - b^{\prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta)\right) = w_i b^{\prime \prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta) (y_i^{\ast}(\boldsymbol \beta) - \mathbf{x}_i^{\prime} \boldsymbol \beta) = \phi w_i(\boldsymbol \beta) (y_i^{\ast}(\boldsymbol \beta) - \mathbf{x}_i^{\prime} \boldsymbol \beta),
\]
a partir de la iteración de Newton-Raphson en la ecuación (13.15), tenemos
\(\boldsymbol \beta_{NEW}\)
\[\begin{eqnarray*}
~ &=& \boldsymbol \beta_{OLD} - \left( \sum_{i=1}^n w_i b^{\prime \prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta _{OLD}) \mathbf{x}_i \mathbf{x}_i^{\prime} \right)^{-1} \left( \sum_{i=1}^n w_i (y_i - b^{\prime}(\mathbf{x}_i^{\prime} \boldsymbol \beta_{OLD})) \mathbf{x}_i \right) \\
&=& \boldsymbol \beta_{OLD} - \left( \sum_{i=1}^n \phi w_i(\boldsymbol \beta _{OLD}) \mathbf{x}_i \mathbf{x}_i^{\prime} \right)^{-1} \left( \sum_{i=1}^n \phi w_i(\boldsymbol \beta _{OLD}) (y_i^{\ast}(\boldsymbol \beta_{OLD}) - \mathbf{x}_i^{\prime} \boldsymbol \beta_{OLD})\mathbf{x}_i \right) \\
&=& \boldsymbol \beta_{OLD} - \left( \sum_{i=1}^n w_i(\boldsymbol \beta _{OLD}) \mathbf{x}_i \mathbf{x}_i^{\prime} \right)^{-1} \left( \sum_{i=1}^n w_i(\boldsymbol \beta _{OLD})\mathbf{x}_i y_i^{\ast}(\boldsymbol \beta_{OLD}) - \sum_{i=1}^n w_i(\boldsymbol \beta _{OLD})\mathbf{x}_i \mathbf{x}_i^{\prime} \boldsymbol \beta_{OLD}) \right) \\
&=& \left( \sum_{i=1}^n w_i(\boldsymbol \beta _{OLD}) \mathbf{x}_i \mathbf{x}_i^{\prime} \right)^{-1} \left( \sum_{i=1}^n w_i(\boldsymbol \beta _{OLD})\mathbf{x}_i y_i^{\ast}(\boldsymbol \beta_{OLD}) \right).
\end{eqnarray*}\]
Así, esto proporciona un método para la iteración usando mínimos cuadrados ponderados. Los mínimos cuadrados reponderados iterativos también están disponibles para enlaces generales utilizando la puntuación de Fisher. Ver McCullagh y Nelder (1989) para más detalles.
