Examen 2

Ejercicio 1 Preferencia de consumidores

Use un modelo linear generalizado Bernoulli/Binomial para analizar la base datos “prefauto.txt” en donde columna 1- americano (1), japonés (0), columna 2 - edad (en años), columna 3 - sexo: masculino (0), femenino (1), columna 4: - estado civil casado (0), soltero (1). Realizar introducción de los datos, análisis descriptivo, inferencial, análisis residual y conclusiones.

Introducción

En este estudio, se busca analizar la preferencia de los consumidores respecto al tipo de automóvil, ya sea americano o japonés, utilizando un modelo lineal generalizado con distribución Bernoulli/Binomial. Se ha tomado una muestra aleatoria de 263 consumidores, y se han registrado varias características sociodemográficas, tales como la edad, el sexo y el estado civil. Las variables observadas incluyen:

• Preferencia por el tipo de automóvil: americano (1) o japonés (0).

• Edad: en años.

• Sexo: masculino (0) o femenino (1).

• Estado civil: casado (0) o soltero (1).

El objetivo es determinar qué factores influyen en la preferencia de los consumidores por automóviles americanos frente a los japoneses, lo que permitirá identificar patrones de comportamiento y posibles correlaciones entre las variables sociodemográficas y las elecciones de compra de automóviles. Para lograr esto, se utilizará un modelo Bernoulli/Binomial que permitirá estimar la probabilidad de que un consumidor prefiera un tipo de automóvil en función de las variables independientes observadas. Este análisis será útil para comprender las tendencias de consumo en este mercado y podría aportar información valiosa para el diseño de estrategias de marketing en la industria automotriz.

De la Tabla 1, Se observa que, independientemente del sexo, la preferencia por los automóviles japoneses es mayoritaria. Entre los individuos casados, aunque la diferencia es leve, también prevalece una inclinación hacia los automóviles japoneses. Sin embargo, esta preferencia es considerablemente más pronunciada entre los solteros, lo que indica una tendencia más clara en este grupo demográfico.

Análisis Descriptivo

prefauto <- read.table("C:/Users/Javier Escobar/Documents/PERSONAL/UNAL/PREGRADO/Asignaturas/2024_1/Modelos Lineales G/Examen2MLG/prefauto.txt", header = FALSE)
# Asignar nombres a las columnas
colnames(prefauto) <- c("Preferencia", "Edad", "Sexo", "EstadoCivil")
head(prefauto)

  Preferencia Edad Sexo EstadoCivil
1           1   34    0           0
2           0   36    0           1
3           0   23    0           0
4           1   29    0           1
5           1   39    0           0
6           0   34    0           1

library(ggplot2)

prefauto$preferencia <- factor(prefauto$Preferencia, levels = c(0, 1), labels = c("Japones", "Americano"))

# Crear el boxplot
ggplot(prefauto, aes(x = preferencia, y = Edad, fill = preferencia)) +
  geom_boxplot() +
  labs(title = "Comparacion de la Edad de los Compradores por Tipo de Automovil",
       x = "Tipo de Automovil",
       y = "Edad") +
  theme_minimal() +
  theme(legend.position = "none") +
  scale_fill_manual(values = c("Japones" = "#9999FF", "Americano" = "#FF9999")) # Personaliza los colores

En los boxplots presentados, se observa que la mediana de edad de los compradores de automóviles americanos es ligeramente superior a la mediana de edad de los compradores de automóviles japoneses. Este hallazgo sugiere que los consumidores que prefieren automóviles americanos tienden a ser un poco mayores en comparación con aquellos que optan por vehículos japoneses. Además, es importante notar que la variabilidad en la edad de los compradores también puede diferir entre los dos grupos, lo que podría indicar distintas características demográficas o preferencias de consumo asociadas a cada tipo de automóvil. Entonces podemos decir que hay una mayor variabilidad en la edad de las personas que prefieren autos americanos (se presenta un dato atipico) es mayor que los que prefieren autos japoneses (se presenta dos dato atipico), y parece que a mayor edad se tiene mayor preferencia a los autos americanos.

Análisis Inferencial

Para analizar la preferencia de los compradores en relación al tipo de automóvil (americano o japonés), es fundamental seleccionar un enfoque estadístico que se ajuste adecuadamente a la naturaleza de los datos. En este contexto, dado que la variable de interés \(Y_i\) es dicotómica (donde \(Y_i = 1\) indica preferencia por automóviles americanos y \(Y_i = 0\) por automóviles japoneses), el modelo de regresión logística se presenta como una opción adecuada.

El modelo logístico permite modelar la probabilidad de un evento binario en función de una o más variables predictoras. Al considerar que \(Y_i\) sigue una distribución Bernoulli, podemos expresar esta relación a través del modelo:

\[ Y_i \sim \text{Bernoulli}(\pi_i), \]

donde \(\pi_i\) representa la probabilidad de que el \(i\)-ésimo comprador prefiera un automóvil americano. Este enfoque no solo facilita la interpretación de los coeficientes estimados, sino que también permite la inclusión de variables explicativas como la edad, el sexo y el estado civil, ofreciendo una comprensión más integral de los factores que influyen en la preferencia de los compradores.

Denotemos por \(Y_i\) la preferencia del \(i\)-ésimo comprador en relación al tipo de automóvil, donde \(1\) indica preferencia por automóviles americanos y \(0\) por automóviles japoneses. Supongamos un modelo logístico sin interacciones, en el que \(Y_i \sim \text{Ber}(\pi_i)\), con la siguiente expresión para el logit de la probabilidad:

\[ \log \left( \frac{\pi_i}{1 - \pi_i} \right) = \beta_0 + \beta_1 \cdot \text{Edad}_i + \beta_2 \cdot \text{Sexo}_i + \beta_3 \cdot \text{EstadoCivil}, \]

donde \(\pi_i\) representa la probabilidad de que el \(i\)-ésimo comprador prefiera un automóvil americano. Este modelo permite evaluar cómo la edad, el sexo y el estado civil influyen en la preferencia por el tipo de automóvil.

Sabemos que \(\pi_i\) es la probabilidad de que el evento ocurra y viene dada por la función logística, la cual es:

\[ \pi_i = \frac{\exp\left\{ \beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2 + x_{3i}\beta_3\right\}}{1 + \exp\left\{ \beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2 + x_{3i}\beta_3 \right\}} \]

A continuación se presenta los resultados de este modelo

Modelo 1

modelo_bernoulli_1 <- glm(preferencia ~ Edad + Sexo + EstadoCivil,
family = binomial(link = "logit"), data = prefauto)
summary(modelo_bernoulli_1)

Call:
glm(formula = preferencia ~ Edad + Sexo + EstadoCivil, family = binomial(link = "logit"), 
    data = prefauto)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -1.65284    0.70771  -2.335   0.0195 *
Edad         0.04982    0.02158   2.309   0.0209 *
Sexo         0.09418    0.25557   0.369   0.7125  
EstadoCivil -0.51787    0.27253  -1.900   0.0574 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 360.44  on 262  degrees of freedom
Residual deviance: 349.70  on 259  degrees of freedom
AIC: 357.7

Number of Fisher Scoring iterations: 4

De lo anterior podemos resaltar que, Las variables Edad y EstadoCivil son las únicas que tienen un impacto estadísticamente significativo en la preferencia de los compradores por automóviles. Entonces por lo que se pudo evidenciar con los p valores, el modelo sugiere que la edad y el estado civil son factores relevantes para la preferencia por automóviles americanos en comparación con los japoneses, mientras que el sexo no parece tener un efecto significativo. por tanto se propone quitar la covariable Sexo del modelo y aplicando el método del criterio de información de Akaike (AIC) para comparar la calidad relativa de los modelos.

Por lo tanto nuestro modelo se especifica como:

\[ \log \left( \frac{\pi_i}{1 - \pi_i} \right) = \beta_0 + \beta_1 \cdot \text{Edad}_i + \beta_2 \cdot \text{EstadoCivil}, \]

Sabemos que \(\pi_i\) es la probabilidad de que el evento ocurra y viene dada por la función logística, la cual es:

\[ \pi_i = \frac{\exp\left\{ \beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2 \right\}}{1 + \exp\left\{ \beta_0 + x_{1i}\beta_1 + x_{2i}\beta_2 \right\}} \]

Modelo 2

modelo_bernoulli_2 <- glm(preferencia ~ Edad + EstadoCivil,
family = binomial(link = "logit"), data = prefauto)
summary(modelo_bernoulli_2)

Call:
glm(formula = preferencia ~ Edad + EstadoCivil, family = binomial(link = "logit"), 
    data = prefauto)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -1.60033    0.69181  -2.313   0.0207 *
Edad         0.04959    0.02153   2.303   0.0213 *
EstadoCivil -0.52578    0.27161  -1.936   0.0529 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 360.44  on 262  degrees of freedom
Residual deviance: 349.83  on 260  degrees of freedom
AIC: 355.83

Number of Fisher Scoring iterations: 4

Antes de proceder a la interpretación de los resultados, es relevante destacar la comparación de los criterios de información de Akaike (AIC) entre los modelos evaluados. Para el modelo 1, el AIC se reportó en 357.7, mientras que para el modelo 2, el AIC fue de 355.83. Esta diferencia indica que el modelo 2 proporciona un mejor ajuste a los datos, ya que un valor de AIC más bajo sugiere un modelo que equilibra mejor la bondad de ajuste y la complejidad. Por lo tanto, podemos concluir que el modelo 2 es preferible para describir la relación entre las variables demográficas y la preferencia por automóviles, lo que justifica su uso en el análisis posterior.

El intercepto ( \(\widehat\beta_0= -1.60033\)) representa el logaritmo de las probabilidades (log-odds) de preferir un automóvil americano cuando todas las demás variables (Edad y Estado Civil) son iguales a cero. Aunque su interpretación directa puede no ser relevante en el contexto, sirve como punto de referencia para los efectos de las otras variables.

El coeficiente para la covariable Edad (\(\widehat\beta_1= 0.04959\)) es positivo indica que, a medida que aumenta la edad de los compradores, también aumenta la probabilidad de que prefieran un automóvil americano. Por cada año adicional de edad, la probabilidad de preferir un automóvil americano en lugar de uno japonés aumenta, dado que el coeficiente es significativo (p = 0.0213). Este efecto sugiere que los compradores más mayores tienden a preferir automóviles americanos.

Para el coeficiente de Estado Civil (\(\widehat\beta_2= -0.52578\)) es negativo sugiere que los compradores casados tienen menor probabilidad de preferir un automóvil americano en comparación con los solteros. Esto sugiere una tendencia que podría ser relevante en el análisis. Ya que implica que el estado civil puede influir en la preferencia por el tipo de automóvil.

Podemos revisar que preferencia tiene mayores probabilidades usando la estimativa de \(\widehat{\pi}\) Así, la probabilidad ajustada de preferencia por automóvil americano queda expresada de la siguiente forma reemplazando las estimaciones de los parámetros de interés descritos anteriormente.

\[ \widehat{\pi} = \frac{\exp(-1.600 + 0.050 \times \text{Edad} - 0.526 \times \text{EstadoCivil})}{1 + \exp(-1.600 + 0.050 \times \text{Edad} - 0.526 \times \text{EstadoCivil})} \]

En el siguiente grafico: Utilizando las estimaciones de \(\widehat{\pi}\), podemos visualizar cómo se comportan las probabilidades de preferencia por automóviles americanos en función de la covariable Estado Civil (Soltero o Casado) a través de un gráfico. Esta representación nos permitirá observar las diferencias en las preferencias de los compradores según su estado civil, facilitando así una comprensión más clara de cómo esta variable influye en la decisión de compra.

library(ggplot2)
library(dplyr)

edades <- seq(min(prefauto$Edad), max(prefauto$Edad), length.out = 100)

# Calcular las probabilidades ajustadas para solteros y casados
probabilidades <- data.frame(
  Edad = rep(edades, 2),
  EstadoCivil = rep(c( "Casado","Soltero"), each = length(edades)),
  Preferencia = c(
    exp(-1.60033 + 0.04959 * edades - 0.52578 * 0) / (1 + exp(-1.60033 + 0.04959 * edades - 0.52578 * 0)), # Casado
    exp(-1.60033 + 0.04959 * edades - 0.52578 * 1) / (1 + exp(-1.60033 + 0.04959 * edades - 0.52578 * 1))  # Soltero
  )
)

# Graficar las probabilidades ajustadas
ggplot(probabilidades, aes(x = Edad, y = Preferencia, color = EstadoCivil)) +
  geom_line(size = 1) +
  labs(title = "Probabilidades Ajustadas de Preferencia por Automoviles Americanos",
       x = "Edad",
       y = "Probabilidad de Preferencia") +
  theme_minimal()

se observa que la preferencia por automóviles americanos aumenta con la edad del comprador. En relación con el estado civil, se nota que los casados prefieren más el automóvil americano que los solteros.

La razon de probalilidades se expresa como

\[ \text{Odds} = \frac{\pi}{1 - \pi} \]

La razón de probabilidades se utiliza para comparar las probabilidades de un evento (en este caso, preferir un automóvil americano) entre dos grupos (casados y solteros).Este mide cuántas veces más probable es que ocurra un evento en un grupo comparado con otro. Se calcula como:

\[ \text{Odds}= \frac{P(\text{Americano} \mid \text{Casado})}{P(\text{Americano} \mid \text{Soltero})}= \exp(\widehat\beta_2)= \exp(0.52578)=1.691778 \]

Esto significa que un comprador casado tiene aproximadamente 1.69 veces más probabilidades de preferir un automóvil americano en comparación con un comprador soltero. Como la razón de probabilidades es mayor que 1, significa que las probabilidades aumentan.

Mientras que una estimación aproximada del intervalo del \(90\%\) para la razón de probabilidades está dada a continuación, y sabemos que el erro estandar para la estimacion de \(\widehat\beta_2\) es \(0.27161\approx 0.272\). \[ e^{0.52578 \pm 1,65 \times 0.272} = e^{0.52578 \pm 0,449} = [1.080; 2.651]. \]

Para interpretar este intervalo en términos de un porcentaje de aumento tenemos que:

Si la razón de probabilidades es mayor que 1, el aumento en probabilidades se calcula como:

\[ (\text{valor} - 1) \times 100 \]

Para el límite superior del intervalo de confianza (2.651), tenemos:

\[ (2.651 - 1) \times 100 = 1.651 \times 100 = 165.1\% \]

Para el límite inferior del intervalo de confianza (1.080), tenemos:

\[ (1.080 - 1) \times 100 = 0.080 \times 100 = 8.0\% \]

Por lo tanto, un comprador casado tiene entre un 8.0% y un 165.1% más de probabilidades de preferir un automóvil americano en comparación con un comprador soltero.

En general las variables Edad y Estado Civil tienen efectos significativos en la preferencia por automóviles americanos. Esto indica que, en el contexto de este análisis, estas características demográficas son importantes para entender las preferencias de los consumidores. por otro parte, dado que Sexo no se incluye en este análisis y las variables restantes muestran significancia, se sugiere mantener solo Edad y Estado Civil en el modelo. La eliminación de variables no significativas nos permite simplificar el modelo y mejorar su interpretación.

Análisis residual

El residuo componente desvío para respuestas Bernoulli está dado por:

\[ TDi = - \frac{2|\ln(1 - \hat{\mu}_i)|^{1/2}}{\sqrt{1 - \hat{h}_{ii}}} \mathbb{I}\{Y_i = 0\} + \frac{2|\ln(\hat{\mu}_i)|^{1/2}}{\sqrt{1 - \hat{h}_{ii}}} \mathbb{I}\{Y_i = 1\} \]

# Programa extraído do site: https://www.ime.usp.br/~giapaula/textoregressao.htm
# Créditos: Prof. Dr. Gilberto Alvarenga Paula
# Adaptado por Caio L. N. Azevedo
# source("D:\\windows\\Unicamp\\Disciplinas\\1_semestre_2016\\ME 720 MLG\\Programas\\diag_Bern.r")

diagBern<-function(fit.model){
  # fit.model: objeto com o resultado do ajuste do MLG obtido através da função glm
  X <- model.matrix(fit.model)
  n <- nrow(X)
  p <- ncol(X)
  w <- fit.model$weights
  W <- diag(w)
  H <- solve(t(X)%*%W%*%X)
  H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
  h <- diag(H)
  ts <- resid(fit.model,type="pearson")/sqrt(1-h)
  td <- resid(fit.model,type="deviance")/sqrt(1-h)
  di <- (h/(1-h))*(ts^2)
  a <- max(td)
  b <- min(td)
  par(mfrow=c(2,2))
  plot(td,xlab="Indice", ylab="Residuo Componente do Desvio",
       ylim=c(b-1,a+1), pch=16,cex.axis=1.1,cex.lab=1.1,cex=1.1,cex.axis=1.1,cex.lab=1.1)
  abline(2,0,lty=2)
  abline(-2,0,lty=2)
  abline(0,0,lty=2)
  
  # identify(td, n=1)
  # title(sub="(c)")
  fited = fitted(fit.model)
  plot(fited ,td,xlab="valor ajustado (media)", ylab="Residuo Componente do Desvio",ylim=c(b-1,a+1), pch=16,
       main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  abline(2,0,lty=2)
  abline(-2,0,lty=2)
  abline(0,0,lty=2)
  
  #
  hist(td,xlab="Residuo Componente do Desvio",ylab="densidad",probability=TRUE,main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  #
  eta = predict(fit.model)
  z = eta + resid(fit.model, type="pearson")/sqrt(w)
  plot(predict(fit.model),z,xlab="Preditor Linear",ylab="Variavel z", pch=16,main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  lines(smooth.spline(predict(fit.model), z, df=2))
  #
  #---------------------------------------------------------------#
}

diagBern(modelo_bernoulli_2)

Del primer gráfico de acuerdo al residuo componente devio (superior izquierda), podemos decir que existe homocedasticidad (la varianza de los errores es constante), ya que los puntos entre los grupos presentan un patron. Además, no advierte de puntos atípicos.

Por otro lado, el gráfico en la parte superior derecha confirma que no se identifica los puntos que pueden ser considerados como atípicos (están por fuera de las líneas de -2 y 2) o mediciones erróneas. Además, podemos ver que los residuos están distribuidos de alrededor de cero y muestran un patrón sistemático.

El histograma de los residuos muestra la distribución de los mismos. En un buen modelo, se espera que los residuos sigan una distribución simétrica alrededor de cero. En nuestro caso, el histograma no es simétrico. debido a la difrencia entre las preferencia de los grupos.

El predictor lineal parece no ajustar bien a los puntos. Esto no respalda la idoneidad del uso del modelo logistico para el análisis de preferencias de automoviles en en función de las variables predictoras.

# Programa extraído do site: https://www.ime.usp.br/~giapaula/textoregressao.htm
# Créditos: Prof. Dr. Gilberto Alvarenga Paula
# Adaptado por Caio L. N. Azevedo
# source("D:\\windows\\Unicamp\\Disciplinas\\1_semestre_2016\\ME 720 MLG\\Programas\\envel_Bern.r")


envelBern <- function(fit.model,ligacao){
  # fit.model: objeto com o resultado do ajuste do MLG obtido através da função glm
  # ligacao: função de ligação (mesmo nome usado pela função glm (colocar entre aspas)
  par(mfrow=c(1,1))
  X <- model.matrix(fit.model)
  n <- nrow(X)
  p <- ncol(X)
  w <- fit.model$weights
  W <- diag(w)
  H <- solve(t(X)%*%W%*%X)
  H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
  h <- diag(H)
  td <- resid(fit.model,type="deviance")/sqrt(1-h)
  e <- matrix(0,n,100)
  #
  for(i in 1:100){
    dif <- runif(n) - fitted(fit.model)
    dif[dif >= 0 ] <- 0
    dif[dif<0] <- 1
    nresp <- dif
    fit <- glm(nresp ~ X, family=binomial(link=ligacao))
    w <- fit$weights
    W <- diag(w)
    H <- solve(t(X)%*%W%*%X)
    H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
    h <- diag(H)
    e[,i] <- sort(resid(fit,type="deviance")/sqrt(1-h))}
  #
  e1 <- numeric(n)
  e2 <- numeric(n)
  #
  for(i in 1:n){
    eo <- sort(e[i,])
    e1[i] <- (eo[2]+eo[3])/2
    e2[i] <- (eo[97]+eo[98])/2}
  #
  med <- apply(e,1,mean)
  faixa <- range(td,e1,e2)
  #par(pty="s")
  qqnorm(td,xlab="Percentil da N(0,1)",
         ylab="Residuo Componente do Desvio", ylim=faixa, pch=16,main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  #
  par(new=T)
  #
  qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1,main="")
  par(new=T)
  qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1,main="")
  par(new=T)
  qqnorm(med,axes=F,xlab="", ylab="", type="l",ylim=faixa,lty=2,main="")
  #---------------------------------------------------------------#                      
}

envelBern(modelo_bernoulli_2, "logit")

Este gráfico refleja lo mencionado anteriormente, por tal razón, vemos que parece haber una tendencia general de los puntos a seguir la línea central, lo que sugiere que los residuos se aproximan a una distribución normal. Pero, se observa que los puntos estan sobre la linea desde -3 a 0 y de 0 a 3, pero en cero ocurre una discontinuidad, es decir no se aproximan ningun punto en esa parte. Esto confirma puede estar ocurriendo problemas con el modelo, aunque la distribución asumida de la variable respuesta es adecuada, quizas no se esta cumpliendo la independencia o la función de enlace no es la adecuada.

Ejercicio 2

Use un modelo linear generalizado Poisson para analizar la base datos “store.txt” en donde columnas : 1 - número de clientes ; 2 - número de hogares ; 3 - renta ; 4 - edad ; 5 - distancia al cliente ; 6 - distancia a la tienda. Realizar introducción de los datos, análisis descriptivo, inferencial, análisis residual y conclusiones.

Introducción

La base de datos “store.txt” contiene seis variables y 110 observaciones sobre el perfil de los clientes de una tienda en 110 áreas de una ciudad. El propósito de este estudio es modelar el número esperado de clientes en cada área utilizando varias variables explicativas como: El número de hogares (en miles), renta anual promedio (en miles de USD), edad promedio de los hogares (en años), distancia al cliente (de otra tienda más cercana que a la tienda en estudio en millas) y distancia a la tienda (en millas). La unidad experimental en este análisis es el área geográfica. Estas variables se recogen con el fin de entender qué factores están asociados con la afluencia de clientes a la tienda.

Dado que la variable dependiente es el número de clientes es un conteo no negativo, el uso de un Modelo Lineal Generalizado (GLM) con distribución Poisson es adecuado para capturar la relación entre las variables predictoras y el número de clientes. Este enfoque permite modelar adecuadamente la naturaleza discreta de los datos y explorar las relaciones entre las variables de interés y la demanda de los clientes en diferentes áreas.

library(ggplot2)
library(dplyr)

store <- read.table("C:\\Users\\Javier Escobar\\Documents\\PERSONAL\\UNAL\\PREGRADO\\Asignaturas\\2024_1\\Modelos Lineales G\\Examen2MLG\\store.txt", header = FALSE)

# Renombrar las columnas para facilitar el análisis
colnames(store) <- c("num_clientes", "num_hogares", "renta", "edad", "distancia_cliente", "distancia_tienda")
# La dimension de la data
dim(store)

[1] 110   6

head(store)

  num_clientes num_hogares renta edad distancia_cliente distancia_tienda
1            9         606 41393    3              3.04             6.32
2            6         641 23635   18              1.95             8.89
3           28         505 55475   27              6.54             2.05
4           11         866 64646   31              1.67             5.81
5            4         599 31972    7              0.72             8.11
6            4         520 41755   23              2.24             6.81

Análisis Descriptivo

# Comprobar si hay valores faltantes
any(is.na(store))

[1] FALSE

# Visualizar la estructura de los datos
str(store)

'data.frame':   110 obs. of  6 variables:
 $ num_clientes     : int  9 6 28 11 4 4 0 14 16 13 ...
 $ num_hogares      : int  606 641 505 866 599 520 354 483 1034 456 ...
 $ renta            : int  41393 23635 55475 64646 31972 41755 46014 34626 85207 33021 ...
 $ edad             : int  3 18 27 31 7 23 26 1 13 32 ...
 $ distancia_cliente: num  3.04 1.95 6.54 1.67 0.72 2.24 0.77 3.51 4.23 3.07 ...
 $ distancia_tienda : num  6.32 8.89 2.05 5.81 8.11 6.81 9.27 7.92 4.4 6.03 ...

# Mostrar un resumen general de los datos
summary(store)

  num_clientes   num_hogares         renta             edad      
 Min.   : 0.0   Min.   :  19.0   Min.   : 19673   Min.   : 1.00  
 1st Qu.: 7.0   1st Qu.: 472.2   1st Qu.: 35160   1st Qu.:13.00  
 Median :10.0   Median : 647.0   Median : 44565   Median :27.00  
 Mean   :11.2   Mean   : 647.8   Mean   : 48837   Mean   :27.43  
 3rd Qu.:14.0   3rd Qu.: 825.2   3rd Qu.: 58369   3rd Qu.:41.75  
 Max.   :32.0   Max.   :1289.0   Max.   :120065   Max.   :58.00  
 distancia_cliente distancia_tienda
 Min.   :0.340     Min.   :0.870   
 1st Qu.:1.927     1st Qu.:5.593   
 Median :2.930     Median :7.280   
 Mean   :3.068     Mean   :6.832   
 3rd Qu.:4.275     3rd Qu.:8.668   
 Max.   :6.610     Max.   :9.900   

El número de clientes varía de 0 a 32 dependiendo del área, y en promedio, hay alrededor de 11 clientes por área.

Hay bastante variación en el número de hogares por área, lo que influirá directamente en el número de clientes potenciales.

La renta promedio varía considerablemente entre las áreas cercanas a la tienda. Las áreas con mayores rentas podrían tener diferentes patrones de compra y visita.

Las áreas alrededor de la tienda tienen una edad promedio que va desde hogares jóvenes hasta mayores, con un promedio de 27 años. Este factor puede influir en los tipos de productos demandados.

La distancia promedio de los clientes a la tienda varía entre áreas, pero la mayoría de los clientes viven a menos de 3 millas. Eso sugiere que la proximidad a la tienda podría ser un factor importante en el número de clientes.

Aunque los hogares están relativamente cerca de la tienda, la distancia puede influir en la decisión de visitar o no la tienda.

# Calcular estadísticas descriptivas básicas para la variable dependiente
descriptivas <- store %>%
  summarise(
    media_clientes = mean(num_clientes),
    var_clientes = var(num_clientes),
    sd_clientes = sd(num_clientes),
  )

# Mostrar las estadísticas calculadas
print(descriptivas)

  media_clientes var_clientes sd_clientes
1           11.2     44.08807    6.639885

El número de clientes varía considerablemente entre las diferentes áreas alrededor de la tienda, con un promedio de 11.2 clientes y una dispersión importante reflejada en una desviación estándar de 6.64. Esto sugiere que hay ciertos factores (como renta, distancia, edad, etc.) que podrían estar afectando el número de clientes de manera desigual entre las áreas estudiadas.

Si suponemos que el número de clientes se distribuye según una distribución de Poisson, entonces la media y la varianza deberían ser aproximadamente iguales, ya que una característica fundamental de la distribución de Poisson es que: \(\text{Var}(Y)=E(Y)\), es decir su varianza y media son iguales.

La media es 11.2, pero la varianza es 44.08, lo que sugiere que los datos no siguen exactamente una distribución de Poisson. En una distribución de Poisson, deberíamos esperar que la varianza sea aproximadamente igual a la media. En este caso, la varianza es considerablemente mayor que la media, lo que indica sobredispersión.

La Sobredispersión ocurre cuando la varianza es mayor que la media. En este caso, sugiere que el modelo de Poisson podría no ser el más adecuado porque subestima la variabilidad observada en los datos. Esto podría deberse a factores no modelados o a una heterogeneidad en los datos.

ggplot(store, aes(x = num_clientes)) +
  geom_histogram(binwidth = 5, fill = "cyan", color = "black") +
  labs(title = "Distribucion del Numero de Clientes", x = "Numero de Clientes", y = "Frecuencia")

La distribución está sesgada a la derecha (positivamente). Esto significa que la mayoría de las áreas tienen un número bajo o moderado de clientes, pero hay algunas áreas que tienen un número considerablemente mayor de clientes (más de 20). Hay una concentración de frecuencias entre 7 y 15 clientes por área, lo que parece ser un rango típico.

Este histograma respalda la hipótesis de sobredispersión en los datos, ya que la distribución no es perfectamente simétrica ni ajusta a una distribución de Poisson simple (donde la varianza debería ser cercana a la media). La presencia de algunas áreas con un número muy elevado de clientes sugiere que podría haber una heterogeneidad en las áreas, con algunas siendo más atractivas para los clientes que otras.

# Crear el boxplot para el número de clientes
boxplot(store$num_clientes, 
        main = "Boxplot del Numero de Clientes", 
        ylab = "Numero de Clientes", 
        col = "cyan", 
        border = "black", 
        horizontal = FALSE)

# Agregar color a la mediana
abline(h = median(store$num_clientes), col = "red")

Hay un posible valor atípico por encima de 25 clientes, indicado por los círculo fuera de los bigotes del boxplot. La distribución parece estar algo sesgada hacia valores mayores, ya que la caja no está completamente centrada.

Veamos la relación del número de clientes con respecto a las demás variables de interes.

# Función para mostrar correlaciones
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) {
    usr <- par("usr"); on.exit(par(usr))  
    par(usr = c(0, 1, 0, 1))  
    r <- round(cor(x, y), digits = digits)  
    txt <- paste0(prefix, r)  
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)  
    text(0.5, 0.5, txt, cex = cex.cor * r)  
}
# Crear la gráfica de pares con correlaciones
pairs(store, 
      col = "green", 
      upper.panel = panel.cor,  # Muestra correlaciones en la parte superior
      lower.panel = points,     # Gráficos de dispersión en la parte inferior
      pch = 16)                 # Tipo de punto en los gráficos

cor(store)

                  num_clientes num_hogares       renta         edad
num_clientes       1.000000000   0.3148981  0.23877229  0.006260648
num_hogares        0.314898135   1.0000000  0.62324747  0.126890363
renta              0.238772295   0.6232475  1.00000000 -0.016555431
edad               0.006260648   0.1268904 -0.01655543  1.000000000
distancia_cliente  0.729556417   0.2572324  0.35720237  0.046250293
distancia_tienda  -0.775625280  -0.2644880 -0.41444964 -0.064488135
                  distancia_cliente distancia_tienda
num_clientes             0.72955642      -0.77562528
num_hogares              0.25723238      -0.26448803
renta                    0.35720237      -0.41444964
edad                     0.04625029      -0.06448814
distancia_cliente        1.00000000      -0.59498958
distancia_tienda        -0.59498958       1.00000000

Observando la primera columna de graficos de disperción se puede concluir que:

1. Número de clientes vs Distancia al cliente: Existe una relación positiva relativamente fuerte (0.73) entre el número de clientes y la distancia del cliente (a otras tiedas competencia de la tienda en estudio). Esto sugiere que, a medida que aumenta la distancia del cliente (en las áreas analizadas), también tiende a aumentar el número de clientes que visitan la tienda. Esta observación es intuitiva, pues normalmente se espera mas clientes visten la tienda dado que las otras tiendas se encuentras un poco más lejos que nuestra tienda de interes.

2. Número de clientes vs Distancia a la tienda: Aquí encontramos una fuerte relación negativa (-0.77), lo cual indica que a medida que aumenta la distancia de la tienda, el número de clientes disminuye. Este comportamiento es esperado, ya que los clientes que viven más lejos de la tienda probablemente tienen menos incentivos para realizar compras allí debido a la lejanía.

3. Número de clientes vs Número de hogares: La correlación positiva (0.31 baja) entre el número de clientes y el número de hogares sugiere que las áreas con más hogares tienden a tener un mayor número de clientes. Aunque la correlación no es muy alta, es lógico esperar que más hogares en una región resulten en más clientes potenciales.

4. Número de clientes vs Renta: La relación entre el número de clientes y la renta es positiva pero débil (0.24 positiva baja), lo que indica que las áreas con mayores rentas tienden a tener más clientes, pero no de manera significativa. Es probable que otros factores además de la renta jueguen un papel más importante en el comportamiento de los clientes.

5. Número de clientes vs Edad: No parece haber ninguna relación significativa entre la edad y el número de clientes (0.006). Esto sugiere que la edad promedio de las áreas no afecta mucho al número de clientes que visitan la tienda.

6. Otras correlaciones destacadas:

Número de hogares vs Renta: Existe una relación positiva moderada (0.62) entre el número de hogares y la renta, lo que podría indicar que las áreas con mayor renta tienden a tener más hogares.

El análisis de correlaciones muestra que las variables más fuertemente relacionadas con el número de clientes son la distancia al cliente (positiva) y la distancia a la tienda (negativa). Estas relaciones sugieren que la proximidad geográfica juega un papel importante en la determinación de la cantidad de clientes, con menos clientes a medida que aumenta la distancia a la tienda y más clientes a medida que aumenta la distancia al cliente de otras tiendas (probablemente debido a la densidad de población o número de hogares).

Análisis inferencial

Se propone el siguiente modelo para nuestro estudio.

Denotamos por \(Y_i\) el número de clientes en el área \(i\)-ésima que fueron a la tienda en un periodo determinado. Suponemos que \(Y_i \sim P(\mu_i)\) con la parte sistemática dada por:

\[ \log(\mu_i) = \beta_0 + \beta_1 \cdot \text{numHogares}_i + \beta_2 \cdot \text{renta}_i + \beta_3 \cdot \text{edad}_i + \beta_4 \cdot \text{distanciaCliente}_i + \beta_5 \cdot \text{distanciaTienda}_i \]

Nota: Cada variable y su respectiva unidad fue definida en la introducción

# Ajustar el modelo Poisson. La variable dependiente es "num_clientes" y las demás son predictores
modelo_poisson1 <- glm(num_clientes ~ num_hogares + renta + edad + distancia_cliente + distancia_tienda,
                      family = poisson(link = "log"), data = store)

# Resumen del modelo
summary(modelo_poisson1)

Call:
glm(formula = num_clientes ~ num_hogares + renta + edad + distancia_cliente + 
    distancia_tienda, family = poisson(link = "log"), data = store)

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)        2.942e+00  2.072e-01  14.198  < 2e-16 ***
num_hogares        6.058e-04  1.421e-04   4.262 2.02e-05 ***
renta             -1.169e-05  2.112e-06  -5.534 3.13e-08 ***
edad              -3.726e-03  1.782e-03  -2.091   0.0365 *  
distancia_cliente  1.684e-01  2.577e-02   6.534 6.39e-11 ***
distancia_tienda  -1.288e-01  1.620e-02  -7.948 1.89e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 422.22  on 109  degrees of freedom
Residual deviance: 114.99  on 104  degrees of freedom
AIC: 571.02

Number of Fisher Scoring iterations: 4

Las variables número de hogares, renta, edad, distancia al cliente y distancia a la tienda son todas estadísticamente significativas (p < 0.05). Esto significa que hay evidencia suficiente para decir que estas variables tienen un efecto en el número de clientes esperado. Los coeficientes estimados en el modelo Poisson con enlace logarítmico representan los efectos de las variables independientes en la media de la variable dependiente, que es el número de clientes.

El valor del intercepto es 2.942, lo que significa que cuando todas las variables predictoras son iguales a 0, el valor esperado del logaritmo del número de clientes es 2.942. Al aplicar la función exponencial (dado que el modelo tiene un enlace logarítmico), el número esperado de clientes sería \(\exp(\widehat{\beta_0})= \exp(2.942)\). Dado que no es realista que estas variables sean todas cero, este valor es principalmente una referencia numérica.

El coeficiente del número de hogares es positivo lo que indica que, a medida que aumenta el número de hogares, también aumenta el número de clientes esperado. Por cada incremento de una unidad en el número de hogares, el logaritmo del número de clientes aumenta en 0.0006058, lo que equivale a un incremento de aproximadamente \(\exp(\widehat{\beta_1})-1= \exp(0.0006058)-1 \thickapprox 0.06\%\) en el número de clientes esperados.

El coeficiente de la renta es negativo, lo que sugiere que a medida que aumenta la renta promedio en un área, el número de clientes esperados disminuye. Por cada incremento de 1000 USD en la renta anual promedio, el número esperado de clientes disminuye en aproximadamente \(\exp(\widehat{\beta_2})-1= \exp(0.00001169)-1 \thickapprox 0.0012\%\).

El coeficiente de la edad también es negativo, lo que sugiere que las áreas con una mayor edad promedio tienen un menor número de clientes esperados. Por cada incremento de un año en la edad promedio, el número de clientes esperados disminuye en aproximadamente \(\exp(\widehat{\beta_3})-1= \exp(−0.003726)-1 \thickapprox 0.37\%\).

En el caso del coeficiente de distancia al cliente (con respecto a otras tiendas distintas a la del estudio) es positivo. Indica que, a medida que aumenta la distancia del cliente, el número de clientes esperados también aumenta. Por cada incremento de una milla en la distancia al cliente, el número de clientes esperados aumenta en \(\exp(\widehat{\beta_4})-1= \exp(0.1684)-1 \thickapprox 18.3\%\).

Por último el coeficiente de la distancia a la tienda es negativo y altamente significativo, lo que indica que a medida que la distancia a la tienda aumenta, el número de clientes esperados disminuye. Por cada incremento de una milla en la distancia a la tienda, el número de clientes esperados disminuye en \(\exp(\widehat{\beta_5})-1= \exp(−0.1288)-1 \thickapprox 12\%\).

Entonces podemos resaltar que: La variable distancia al competidor tiene un efecto positivo importante en el número de clientes, mientras que la distancia a la tienda tiene un efecto negativo significativo.

La devianza del modelo ajustado (Residual deviance: 114.99 con 104 grados de libertad). Indica una reducción sustancial en la devianza sugiere que el modelo ajustado con predictores mejora el ajuste con respecto al modelo nulo. \[ D(\mathbf{y}; \hat{\boldsymbol{\mu}}) = 2 \sum_{i=1}^n \left\{ y_i \log\left(\frac{y_i}{\hat{\mu}_i}\right) - (y_i - \hat{\mu}_i) \right\}. \]

Análisis residual.

# Creditos: Prof. Dr. Gilberto Alvarenga Paula
# Adaptado por Caio L. N. Azevedo

diagPoisson<-function(fit.model){
  # fit.model: objeto com o resultado do ajuste do MLG obtido através da função glm
  X <- model.matrix(fit.model)
  n <- nrow(X)
  p <- ncol(X)
  w <- fit.model$weights
  W <- diag(w)
  H <- solve(t(X)%*%W%*%X)
  H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
  h <- diag(H)
  ts <- resid(fit.model,type="pearson")/sqrt(1-h)
  td <- resid(fit.model,type="deviance")/sqrt(1-h)
  di <- (h/(1-h))*(ts^2)
  par(mfrow=c(2,2))
  a <- max(td)
  b <- min(td)
  plot(td,xlab="Indice", ylab="Residuo Componente do Desvio",
       ylim=c(b-1,a+1), pch=16,cex.axis=1.1,cex.lab=1.1,cex=1.1,cex.axis=1.1,cex.lab=1.1)
  abline(2,0,lty=2)
  abline(-2,0,lty=2)
  abline(0,0,lty=2)
  
  # identify(td, n=1)
  # title(sub="(c)")
  fited = fitted(fit.model)
  plot(fited ,td,xlab="valor ajustado (media)", ylab="Residuo Componente do Desvio",ylim=c(b-1,a+1), pch=16,
       main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  abline(2,0,lty=2)
  abline(-2,0,lty=2)
  abline(0,0,lty=2)
  
  #
  hist(td,xlab="Residuo Componente do Desvio",ylab="densidad",probability=TRUE,main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  #
  eta = predict(fit.model)
  z = eta + resid(fit.model, type="pearson")/sqrt(w)
  plot(predict(fit.model),z,xlab="Preditor Linear",ylab="Variavel z", pch=16,main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  abline(0,1)
  par(mfrow=c(1,1))
  #------------------------------------------------------------#
  return(td)
}


diagPoisson(modelo_poisson1)

           1            2            3            4            5            6 
-1.017747721 -1.008560368 -0.025678465  0.870515499 -1.350777365 -1.854025104 
           7            8            9           10           11           12 
-2.968627052  0.887800967  0.415185057  0.397487892  0.212061278  0.608740939 
          13           14           15           16           17           18 
-0.294635496 -0.394504308  0.290209855 -0.464215636 -1.039612603  0.148434300 
          19           20           21           22           23           24 
 1.636114701  1.917802139  0.703385909  0.525902451 -1.420497402 -0.723867582 
          25           26           27           28           29           30 
 0.089108198 -0.676142578 -1.252882692  0.309440376  2.278914070 -1.266654342 
          31           32           33           34           35           36 
-2.052913483 -1.297883086 -1.095982149 -0.583017964  0.596585029  1.218238643 
          37           38           39           40           41           42 
-2.924883690  1.988231730  0.785557519  0.415742594 -0.213054768 -0.398029616 
          43           44           45           46           47           48 
-1.750735071 -1.641297116 -2.921106519 -0.213803379 -0.181585868  0.445356527 
          49           50           51           52           53           54 
 0.620824909 -0.613089398 -2.397954691 -0.230873325  0.386723369 -0.005976834 
          55           56           57           58           59           60 
 0.319560473  1.251134016 -0.331631772  0.569516680 -1.627005079  1.827006514 
          61           62           63           64           65           66 
-0.365530858  0.496962815  1.000797516  0.743229550 -0.497417361 -0.225436823 
          67           68           69           70           71           72 
-0.079436149 -0.567773737  0.918914001  1.016118106  0.626514138  0.953456017 
          73           74           75           76           77           78 
-0.196759702 -1.397802039  0.538631188 -0.425721843 -0.764650023 -0.598323137 
          79           80           81           82           83           84 
 0.189959085  1.569328206  1.453545806 -0.281069048  0.021105807  0.382825998 
          85           86           87           88           89           90 
-1.705330238  0.093779903  0.537597535  0.036141931  0.006337590 -1.232647913 
          91           92           93           94           95           96 
-0.997584447  0.613922288  0.836989310  0.299612321  0.132462824  0.162762617 
          97           98           99          100          101          102 
-0.862435832 -0.665170643 -0.385727142 -0.318441611  0.686107748  0.817360577 
         103          104          105          106          107          108 
 1.904107310  1.455208044 -0.594201343  0.818232760 -0.246596386  0.293881449 
         109          110 
-0.200897412 -0.173041057 

Del primer gráfico (superior izquierda), podemos decir que existe heterocedasticidad (la varianza de los errores no es constante), ya que los puntos inician más dispersos y, hacia el final, se reduce la dispersión (están distribuidos aleatoriamente). Además, nos advierte de algunos puntos atípicos.

Por otro lado, el gráfico en la parte superior derecha revela exactamente los puntos que pueden ser considerados como atípicos (están por fuera de las líneas de -2 y 2) o mediciones erróneas. Además, podemos ver que los residuos están distribuidos de manera aleatoria alrededor de cero y no muestran un patrón sistemático.

El histograma de los residuos muestra la distribución de los mismos. En un buen modelo, se espera que los residuos sigan una distribución simétrica alrededor de cero. En nuestro caso, el histograma parece relativamente simétrico, lo que sugiere que los residuos se distribuyen adecuadamente en torno a la media. No obstante, algunos valores más extremos hacia la izquierda podrían sugerir la presencia de posibles outliers o problemas de ajuste en ciertas observaciones.

El predictor lineal parece ajustar bien a los puntos. Aunque algunos valores residuales más grandes pueden señalar posibles puntos atípicos o falta de ajuste para ciertas observaciones, en términos generales no se observan problemas significativos. Esto respalda la idoneidad del uso del modelo de regresión Poisson para el análisis del número de clientes en función de las variables predictoras.

Graficas de envolventes

envelPoisson <- function(fit.model,ligacao){
  # fit.model: objeto com o resultado do ajuste do MLG obtido através da função glm
  # ligacao: função de ligação (mesmo nome usado pela função glm (colocar entre aspas)
  #par(mfrow=c(1,1))
  X <- model.matrix(fit.model)
  n <- nrow(X)
  p <- ncol(X)
  w <- fit.model$weights
  W <- diag(w)
  H <- solve(t(X)%*%W%*%X)
  H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
  h <- diag(H)
  td <- resid(fit.model,type="deviance")/sqrt((1-h))
  e <- matrix(0,n,100)
  #
  for(i in 1:100){
    nresp <- rpois(n, fitted(fit.model))
    fit <- glm(nresp ~ X, family=poisson(link=ligacao))
    w <- fit$weights
    W <- diag(w)
    H <- solve(t(X)%*%W%*%X)
    H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
    h <- diag(H)
    e[,i] <- sort(resid(fit,type="deviance")/sqrt(1-h))}
  #
  e1 <- numeric(n)
  e2 <- numeric(n)
  #
  for(i in 1:n){
    eo <- sort(e[i,])
    e1[i] <- (eo[2]+eo[3])/2
    e2[i] <- (eo[97]+eo[98])/2}
  #
  med <- apply(e,1,mean)
  faixa <- range(td,e1,e2)
  #par(pty="s")
  qqnorm(td,xlab="Percentil da N(0,1)",
         ylab="Residuo Componente do Desvio", ylim=faixa, pch=16,main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  #
  par(new=T)
  #
  qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1,main="")
  par(new=T)
  qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1,main="")
  par(new=T)
  qqnorm(med,axes=F,xlab="", ylab="", type="l",ylim=faixa,lty=2,main="")
  #------------------------------------------------------------#                      
}

envelPoisson(modelo_poisson1, "log")

Este gráfico refleja lo mencionado en el análisis descriptivo acerca de la sobredispersión, ya que la media y la varianza de los datos difieren altamente, supuesto que se debe cumplir para suponer que nuestras respuestas provienen de una distribución Poisson. Por tal razón, vemos que parece haber una tendencia general de los puntos a seguir la línea diagonal, lo que sugiere que los residuos se aproximan a una distribución normal. Sin embargo, se observa que los puntos se alejan significativamente de la línea en los extremos (es decir, en las colas)

Ejercicio 3

Resuelva los puntos. 8.2, 9.1 y 9.2 del Dobson, A. J., Barnett, A. G. (2018). An Introduction to Generalized Linear Models. Estados Unidos: CRC Press.

Punto 9.1

Sea \(Y_1, \ldots, Y_N\) variables aleatorias independientes con \(Y_i \sim \text{Poisson}(\mu_i)\) y

\[ \log \mu_i = \beta_1 + \sum_{j=2}^{J} x_{ij} \beta_j, \quad i = 1, \ldots, N. \]

1. Demuestra que la estadística de puntuación para \(\beta_1\) es

\[ U_1 = \sum_{i=1}^{N} (Y_i - \mu_i). \]

La función de masa de probabilidad de una distribución de Poisson está dada por:

\[ f(y_i) = \frac{\mu_{i}^{y_i} e^{-\mu_{i}}}{y_i !} \]

donde:

• \(y_i\) es el valor observado, \(y_i = 0, 1, 2, \ldots\)
• \(\mu_i\) es el parámetro de media (y varianza) de la distribución de Poisson. \(\mu\) pertenece al intervalo \((0, \infty)\)

Como dicha distribución pertenece a la famila exponencial lineal, la función de densidad se puede escribir como: \[ f(y_i)= \exp\{ y_i \log(\mu_i) - \mu_i - \log(y_i !) \} \] La función de log-verosimilitud, sabiendo que tenemos variables aleatorias independientes se puede expresar como:

\[ l(\mu_{i};y_{i})= \log \prod_{i=1}^{N} f(y_i)= \sum_{i=1}^{N} \log f(y_i)\\ = \sum_{i=1}^{N} \left [ y_i \log(\mu_i) - \mu_i - \log(y_i !) \right ] \]

Por la definición al inicio, \(\mu_i =\exp \left [ \beta_1 + \sum_{j=2}^{J} x_{ij} \beta_j \right ]\)

Sabemos que

\[ U_j = \frac{\partial l}{\partial \beta_j} \Rightarrow \frac{\partial l}{\partial \beta_1} = U_1 = \frac{\partial l_i}{\partial \beta_1} = \left[ \frac{\partial l_i}{\partial \mu_i} \cdot \frac{\partial \mu_i}{\partial \beta_1} \right]\\ \text{Calculando las respectivas derivadas obtenemos}\\ \Rightarrow \frac{\partial l}{\partial \mu_i} = \sum_{i=1}^{N} \left [ y_i (1/\mu_i) - 1 - 0 \right ] \\ \Rightarrow \frac{\partial \mu_i}{\partial \beta_1} = \mu_i =\exp \left [ \beta_1 + \sum_{j=2}^{J} x_{ij} \beta_j \right ]\\ \text{Remplazando logramos mostrar que:} \\ \Rightarrow U_1 = \sum_{i=1}^{N} \left [ y_i (1/\mu_i) - 1 \right ] \cdot \mu_i = \sum_{i=1}^{N} \left [ y_i -\mu_i \right ] \]

3. Por lo tanto, muestra que para las estimaciones de máxima verosimilitud \(\hat{\mu}_i\),

\[ \sum \hat{\mu}_i = \sum y_i. \]

Sabemos que para la hallar el estimador de Máxima verosimilitud, basta con derivar e igualar a cero la función de log-verosimilitud para encontrar puntos críticos y evaluar cual de ellos es quien maximiza la función de log-verosimilitud con el criterio de la segunda derivada, como esta funcion fue construida anteriormente, tenemos que:

\[ \frac{\partial l(\mu_{i};y_{i})}{\partial \mu_i} = \sum_{i=1}^{N} \left [ y_i(1/\mu_i) - 1 \right ].\\ \text{ Entonces:}\\ \sum_{i=1}^{N} \left [ y_i (1/ \ \widehat\mu_i) - 1 \right ]= 0\\ \Leftrightarrow \sum_{i=1}^{N} \frac{y_i}{\widehat\mu_i} - \sum_{i=1}^{N}1 = 0\\ \Leftrightarrow \sum_{i=1}^{N} \frac{y_i}{\widehat\mu_i} = \sum_{i=1}^{N}1 \text{ ; multiplicando en ambos lados de la igualdad } \widehat\mu_i.\\ \Rightarrow \sum_{i=1}^{N} y_i = \sum_{i=1}^{N} \widehat\mu_i. \]

4. Deduce que la expresión para la desviación en (9.6) \(D = 2 \sum \left[ o_i \log\left(\frac{o_i}{e_i}\right) - (o_i - e_i) \right]\) se simplifica a (9.7) \(D = 2 \sum \left[ o_i \log\left(\frac{o_i}{e_i}\right) \right]\) en este caso.

Para entender la solución debemos recordar que:

\[ D = 2 \sum_{i=1}^{N} \left[ y_i \log\left(\frac{y_i}{\widehat\mu_i}\right) - (y_i - \widehat\mu_i) \right]\\ \text{Es decir que: } o_i=y_i \text{ y tambien que: } e_i= \widehat{y_i} = \widehat{\mu}_i = e^{\boldsymbol{x_i^T} \boldsymbol{\widehat{\beta}}} \]

Por lo demostrado en el ejercicio anterior, podemos escribir que: ${i=1}^{N} o_i = {i=1}^{N} e_i $

\[ \text{Remplazando las operaciones necesarias podemos simplicar } \textbf{D :}\\ D = 2 \sum_{i=1}^{N} \left[ o_i \log\left(\frac{o_i}{e_i}\right) - (o_i - e_i) \right]\\ = 2 \sum_{i=1}^{N} \left[ o_i \log\left(\frac{o_i}{e_i}\right) \right] - \sum_{i=1}^{N}(o_i - e_i)\\ = 2 \sum_{i=1}^{N} \left[ o_i \log\left(\frac{o_i}{e_i}\right) \right] - \sum_{i=1}^{N} o_i + \sum_{i=1}^{N} e_i \\ = 2 \sum_{i=1}^{N} \left[ o_i \log\left(\frac{o_i}{e_i}\right) \right] \]

Punto 9.2

Los datos en la Tabla 9.13 son números de pólizas de seguro, \(n\), y números de reclamaciones, \(y\), para automóviles en varias categorías de seguros, , tabulados por edad del titular de la póliza, , y el distrito donde vivía el titular de la póliza ( = 1, para Londres y otras grandes ciudades, y = 0, de lo contrario). La tabla se deriva del conjunto de datos en Aitkin et al. (1989), obtenido de un artículo de Baxter, Coutts y Ross (1980).

1. Calcular la tasa de reclamaciones \(y/n\) para cada categoría y graficar las tasas por , y para tener una idea de los principales efectos de estos factores.

library(dobson)
Tabla_9.13 <- as.data.frame(insurance)
dim(insurance)

[1] 32  5

head(Tabla_9.13)

  car age district   y    n
1   1   1        0  65  317
2   1   2        0  65  476
3   1   3        0  52  486
4   1   4        0 310 3259
5   2   1        0  98  486
6   2   2        0 159 1004

# Para calcular la tasa de reclamaciones
insurance$rate <- insurance$y / insurance$n

# Convertir 'car' a factor para asegurar colores discretos
insurance$car <- as.factor(insurance$car)

library(ggplot2)

# Gráfico de tasas por edad (AGE)
ggplot(insurance, aes(x = age, y = rate, color = car)) +
    geom_point() +
    geom_line() +
    labs(title = "Tasa de Reclamaciones por Edad y Categoria de Seguro",
         x = "Edad del Asegurado",
         y = "Tasa de Reclamaciones (y/n)") +
    facet_wrap(~district) +
    scale_color_manual(values = c("red", "blue", "green", "orange")) +  # Colores distintos para cada categoría
    theme_minimal()

De nuestros datos sabemos que tenemos 4 (1:4) categorias de Edad del asegurado y 4 (1:4) categorias de categorias de seguro. Además tenemos dos distritos (1 = Otros, 0 = Londres y Ciudades Principales). Del anterior grafico podemos decir que:

2. Utilizar la regresión Poisson para estimar los principales efectos (cada uno tratado como categórico y modelado usando variables indicadoras) y términos de interacción.

3. Basado en el modelado en (b), Aitkin et al. (1989) determinaron que todas las interacciones eran insignificantes y decidieron que y podrían tratarse como si fueran variables continuas. Ajustar un modelo que incorpore estas características y compararlo con el mejor modelo obtenido en (b). ¿A qué conclusiones llegas?

library(ggplot2)

# Ajustar el modelo de regresión de Poisson, teniendo en cuenta las 
# interacciones entre todas la covavariables
modelo1 <- glm(y ~ factor(car) * factor(age) * factor(district), 
               family = poisson(link = "log"), 
               data = insurance)

# Resumen del modelo
summary(modelo1)

Call:
glm(formula = y ~ factor(car) * factor(age) * factor(district), 
    family = poisson(link = "log"), data = insurance)

Coefficients:
                                              Estimate Std. Error z value
(Intercept)                                  4.174e+00  1.240e-01  33.655
factor(car)2                                 4.106e-01  1.600e-01   2.567
factor(car)3                                -4.608e-01  1.994e-01  -2.311
factor(car)4                                -1.776e+00  3.260e-01  -5.449
factor(age)2                                 6.624e-15  1.754e-01   0.000
factor(age)3                                -2.231e-01  1.861e-01  -1.199
factor(age)4                                 1.562e+00  1.364e-01  11.451
factor(district)1                           -3.481e+00  7.179e-01  -4.849
factor(car)2:factor(age)2                    4.839e-01  2.174e-01   2.226
factor(car)3:factor(age)2                    1.049e+00  2.524e-01   4.155
factor(car)4:factor(age)2                    1.157e+00  3.876e-01   2.986
factor(car)2:factor(age)3                    8.030e-01  2.248e-01   3.572
factor(car)3:factor(age)3                    1.430e+00  2.575e-01   5.552
factor(car)4:factor(age)3                    1.489e+00  3.888e-01   3.829
factor(car)2:factor(age)4                    6.294e-01  1.731e-01   3.636
factor(car)3:factor(age)4                    8.918e-01  2.124e-01   4.199
factor(car)4:factor(age)4                    1.158e+00  3.399e-01   3.407
factor(car)2:factor(district)1               8.422e-01  8.176e-01   1.030
factor(car)3:factor(district)1               1.377e+00  8.601e-01   1.601
factor(car)4:factor(district)1              -2.122e+01  4.225e+04  -0.001
factor(age)2:factor(district)1               9.163e-01  8.549e-01   1.072
factor(age)3:factor(district)1               9.163e-01  8.858e-01   1.034
factor(age)4:factor(district)1               1.328e+00  7.392e-01   1.797
factor(car)2:factor(age)2:factor(district)1 -1.044e+00  9.950e-01  -1.049
factor(car)3:factor(age)2:factor(district)1 -1.628e+00  1.052e+00  -1.548
factor(car)4:factor(age)2:factor(district)1  2.202e+01  4.225e+04   0.001
factor(car)2:factor(age)3:factor(district)1 -3.510e-01  9.944e-01  -0.353
factor(car)3:factor(age)3:factor(district)1 -9.595e-01  1.039e+00  -0.924
factor(car)4:factor(age)3:factor(district)1  2.220e+01  4.225e+04   0.001
factor(car)2:factor(age)4:factor(district)1 -8.407e-01  8.428e-01  -0.997
factor(car)3:factor(age)4:factor(district)1 -1.248e+00  8.881e-01  -1.406
factor(car)4:factor(age)4:factor(district)1  2.175e+01  4.225e+04   0.001
                                            Pr(>|z|)    
(Intercept)                                  < 2e-16 ***
factor(car)2                                0.010267 *  
factor(car)3                                0.020856 *  
factor(car)4                                5.07e-08 ***
factor(age)2                                1.000000    
factor(age)3                                0.230388    
factor(age)4                                 < 2e-16 ***
factor(district)1                           1.24e-06 ***
factor(car)2:factor(age)2                   0.026013 *  
factor(car)3:factor(age)2                   3.26e-05 ***
factor(car)4:factor(age)2                   0.002826 ** 
factor(car)2:factor(age)3                   0.000354 ***
factor(car)3:factor(age)3                   2.83e-08 ***
factor(car)4:factor(age)3                   0.000129 ***
factor(car)2:factor(age)4                   0.000277 ***
factor(car)3:factor(age)4                   2.68e-05 ***
factor(car)4:factor(age)4                   0.000657 ***
factor(car)2:factor(district)1              0.302970    
factor(car)3:factor(district)1              0.109356    
factor(car)4:factor(district)1              0.999599    
factor(age)2:factor(district)1              0.283777    
factor(age)3:factor(district)1              0.300931    
factor(age)4:factor(district)1              0.072362 .  
factor(car)2:factor(age)2:factor(district)1 0.294295    
factor(car)3:factor(age)2:factor(district)1 0.121616    
factor(car)4:factor(age)2:factor(district)1 0.999584    
factor(car)2:factor(age)3:factor(district)1 0.724125    
factor(car)3:factor(age)3:factor(district)1 0.355571    
factor(car)4:factor(age)3:factor(district)1 0.999581    
factor(car)2:factor(age)4:factor(district)1 0.318567    
factor(car)3:factor(age)4:factor(district)1 0.159809    
factor(car)4:factor(age)4:factor(district)1 0.999589    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 5.6606e+03  on 31  degrees of freedom
Residual deviance: 4.1217e-10  on  0  degrees of freedom
AIC: 232.36

Number of Fisher Scoring iterations: 20

# Modelo de interacciones
modelo2 <- glm(y ~ car * age * district, 
               family = poisson(link = "log"), 
               data = insurance)
summary(modelo2)

Call:
glm(formula = y ~ car * age * district, family = poisson(link = "log"), 
    data = insurance)

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)        2.91917    0.16066  18.170  < 2e-16 ***
car2               0.30090    0.19688   1.528 0.126416    
car3              -0.05176    0.21492  -0.241 0.809675    
car4              -1.53013    0.31386  -4.875 1.09e-06 ***
age                0.65565    0.04769  13.749  < 2e-16 ***
district          -4.02026    0.75664  -5.313 1.08e-07 ***
car2:age           0.20416    0.05770   3.539 0.000402 ***
car3:age           0.15227    0.06296   2.418 0.015586 *  
car4:age           0.25783    0.08996   2.866 0.004158 ** 
car2:district      0.94067    0.87021   1.081 0.279709    
car3:district      1.14972    0.91019   1.263 0.206529    
car4:district      1.53013    1.09174   1.402 0.161049    
age:district       0.48893    0.20832   2.347 0.018924 *  
car2:age:district -0.24427    0.23979  -1.019 0.308341    
car3:age:district -0.27576    0.25172  -1.095 0.273298    
car4:age:district -0.25783    0.30056  -0.858 0.390989    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 5660.59  on 31  degrees of freedom
Residual deviance:  285.39  on 16  degrees of freedom
AIC: 485.75

Number of Fisher Scoring iterations: 5

# modelo simple 
modelo3 <- glm(y ~ car + age + district, 
               family = poisson(link = "log"), 
               data = insurance)
summary(modelo3)

Call:
glm(formula = y ~ car + age + district, family = poisson(link = "log"), 
    data = insurance)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.30209    0.08232  27.965  < 2e-16 ***
car2         0.98960    0.05045  19.617  < 2e-16 ***
car3         0.47070    0.05490   8.574  < 2e-16 ***
car4        -0.58927    0.07211  -8.172 3.03e-16 ***
age          0.83664    0.02067  40.482  < 2e-16 ***
district    -2.15937    0.05849 -36.917  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 5660.59  on 31  degrees of freedom
Residual deviance:  324.84  on 26  degrees of freedom
AIC: 505.2

Number of Fisher Scoring iterations: 4

#  Modelo simple Teniendo en cuenta el desplazamiento n (offset)
modelo4 <- glm(y ~  car + age + district  +
offset(log(n)),
               family = poisson(link = "log"), 
               data = insurance)
summary(modelo4)

Call:
glm(formula = y ~ car + age + district + offset(log(n)), family = poisson(link = "log"), 
    data = insurance)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.63733    0.07499 -21.833  < 2e-16 ***
car2         0.16260    0.05048   3.221 0.001276 ** 
car3         0.39389    0.05491   7.174 7.31e-13 ***
car4         0.56585    0.07216   7.842 4.44e-15 ***
age         -0.17616    0.01850  -9.523  < 2e-16 ***
district     0.21860    0.05853   3.735 0.000188 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 207.833  on 31  degrees of freedom
Residual deviance:  23.832  on 26  degrees of freedom
AIC: 204.19

Number of Fisher Scoring iterations: 4

Por medio del método del AIC, podemos elegir cual modelo presenta un mejor ajuste a nuestros datos, para este caso vemos que el Modelo4 ” Modelo simple Teniendo en cuenta el desplazamiento n (offset)” presenta el menor AIC en comparacion a los demas modelos, por lo tanto será tomado como el mejor en este caso.

Algunas interpretaciones

El intercepto es significativamente negativo, lo que indica que, cuando todas las variables son cero, la tasa de reclamaciones es baja.

Efectos de las Categorías de Seguro: car2, car3, car4: Todos estos coeficientes son positivos y significativos, lo que indica que a medida que se incrementa la categoría de seguro (comparado con la categoría de referencia), la tasa de reclamaciones aumenta. Esto respalda la conclusión de que las tasas de reclamaciones aumentan con la categoría del seguro.

Edad: El coeficiente de age es negativo y significativo, lo que indica que a medida que aumenta la edad del asegurado, la tasa de reclamaciones disminuye.

Distrito: El coeficiente de district es positivo y significativo, lo que sugiere que los asegurados en distritos de mayor riesgo tienen tasas de reclamaciones más altas.

Diagnóstico del Modelo: El AIC es 204.19, lo que sugiere que es más eficiente que otros modelos considerados.

Desviación Residual: La reducción en la desviación residual también sugiere que tu modelo se ajusta bien a los datos. La comparación de la desviación nula con la desviación residual muestra que el modelo explica una buena parte de la variabilidad en los datos.

El modelo sugiere que las tasas de reclamaciones son influenciadas positivamente por las categorías de seguro y negativamente por la edad del asegurado, además de verse afectadas por el distrito de residencia. Estos hallazgos son consistentes con las conclusiones previas y pueden ser útiles para la evaluación de riesgos en la fijación de precios de pólizas de seguros.

Analisis Residuos Componentes Desvio

diagPoisson<-function(fit.model){
  # fit.model: objeto com o resultado do ajuste do MLG obtido através da função glm
  X <- model.matrix(fit.model)
  n <- nrow(X)
  p <- ncol(X)
  w <- fit.model$weights
  W <- diag(w)
  H <- solve(t(X)%*%W%*%X)
  H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
  h <- diag(H)
  ts <- resid(fit.model,type="pearson")/sqrt(1-h)
  td <- resid(fit.model,type="deviance")/sqrt(1-h)
  di <- (h/(1-h))*(ts^2)
  par(mfrow=c(2,2))
  a <- max(td)
  b <- min(td)
  plot(td,xlab="Indice", ylab="Residuo Componente do Desvio",
       ylim=c(b-1,a+1), pch=16,cex.axis=1.1,cex.lab=1.1,cex=1.1,cex.axis=1.1,cex.lab=1.1)
  abline(2,0,lty=2)
  abline(-2,0,lty=2)
  abline(0,0,lty=2)
  
  # identify(td, n=1)
  # title(sub="(c)")
  fited = fitted(fit.model)
  plot(fited ,td,xlab="valor ajustado (media)", ylab="Residuo Componente do Desvio",ylim=c(b-1,a+1), pch=16,
       main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  abline(2,0,lty=2)
  abline(-2,0,lty=2)
  abline(0,0,lty=2)
  
  #
  hist(td,xlab="Residuo Componente do Desvio",ylab="densidad",probability=TRUE,main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  #
  eta = predict(fit.model)
  z = eta + resid(fit.model, type="pearson")/sqrt(w)
  plot(predict(fit.model),z,xlab="Preditor Linear",ylab="Variavel z", pch=16,main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  abline(0,1)
  par(mfrow=c(1,1))
  #------------------------------------------------------------#
  return(td)
}

# library(glmtoolbox)


diagPoisson(modelo4)

          1           2           3           4           5           6 
 1.97774258 -0.01206370 -0.53376179 -0.31368050  0.56326992 -0.22665661 
          7           8           9          10          11          12 
-0.62657680  0.69674240 -2.01390891  0.81543763  1.79562040 -1.05678768 
         13          14          15          16          17          18 
-0.14948860 -0.11589443 -0.70566313 -0.67070646 -1.14972968 -0.26856083 
         19          20          21          22          23          24 
-0.76576173 -0.32625864 -0.15215775 -1.72239593  0.34478735  0.01216246 
         25          26          27          28          29          30 
-0.18414630 -0.97747503  0.43068413  0.29046933 -1.47072096  0.53863427 
         31          32 
 0.67072327  1.87614125 

En el primer grafico se puede notar la heterocedasticidad se observa un patrón (como un embudo). Ademas parece haber datos atipicos, como se puede notar en la segunda grafica los puntos que sobresalen o estan solo las lineas de (-2, 2). Por tro lado no se observa una simetria clara en la distribución de los residuos, esto tambien advierte de un mal ajuste debido a que en un buen modelo se espera que los residuos sigan una distribución simétrica alrededor de cero. Por ultimo parece que el predictor lineal se ajusta bien a los datos, solo se escapa un punto que podria considerarse atipico.

envelPoisson <- function(fit.model,ligacao){
  # fit.model: objeto com o resultado do ajuste do MLG obtido através da função glm
  # ligacao: função de ligação (mesmo nome usado pela função glm (colocar entre aspas)
  #par(mfrow=c(1,1))
  X <- model.matrix(fit.model)
  n <- nrow(X)
  p <- ncol(X)
  w <- fit.model$weights
  W <- diag(w)
  H <- solve(t(X)%*%W%*%X)
  H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
  h <- diag(H)
  td <- resid(fit.model,type="deviance")/sqrt((1-h))
  e <- matrix(0,n,100)
  #
  for(i in 1:100){
    nresp <- rpois(n, fitted(fit.model))
    fit <- glm(nresp ~ X, family=poisson(link=ligacao))
    w <- fit$weights
    W <- diag(w)
    H <- solve(t(X)%*%W%*%X)
    H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
    h <- diag(H)
    e[,i] <- sort(resid(fit,type="deviance")/sqrt(1-h))}
  #
  e1 <- numeric(n)
  e2 <- numeric(n)
  #
  for(i in 1:n){
    eo <- sort(e[i,])
    e1[i] <- (eo[2]+eo[3])/2
    e2[i] <- (eo[97]+eo[98])/2}
  #
  med <- apply(e,1,mean)
  faixa <- range(td,e1,e2)
  #par(pty="s")
  qqnorm(td,xlab="Percentil da N(0,1)",
         ylab="Residuo Componente do Desvio", ylim=faixa, pch=16,main="",cex=1.1,cex.axis=1.1,cex.lab=1.1)
  #
  par(new=T)
  #
  qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1,main="")
  par(new=T)
  qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1,main="")
  par(new=T)
  qqnorm(med,axes=F,xlab="", ylab="", type="l",ylim=faixa,lty=2,main="")
  #------------------------------------------------------------#                      
}

envelPoisson(modelo4, "log")

Nuestras variable de interes que corresponde al numero de reclamaciones, por ser datos de conteo se asume que viene dada de una distribución Poisson. Entonces la media y la varianza deberían ser aproximadamente iguales, ya que una característica fundamental de la distribución de Poisson es que: \(Var(Y)=E(Y)\), es decir su varianza y media son iguales.

mean(insurance$y)

[1] 98.46875

var(insurance$y)

[1] 30384.13

Como se acaba de mostrar esto no se cumple, lo que sugiere que los datos tienen una mayor variabilidad de lo que sería típico en una distribución de Poisson, que asume que la media y la varianza son iguales.

Esto nos lleva a tener sobredispersión de nuestro modelo. La sobredispersión puede hacer que un modelo de Poisson no se ajuste bien a los datos, lo que podría llevar a estimaciones inexactas de los parámetros y a inferencias incorrectas, porque subestima la variabilidad observada en los datos.

Lo anterior explica lo que esta ocurriendo en la gráfica de envolventes, en donde se puede observar que muy pocos puntos se aproximan a la linea central, en cambio estos sobre salen o estan por fuera de las envolventes.

Tambien se podria considerar el uso de modelos alternativos que manejan mejor la sobredispersión, como el Binomial Negativo este modelo es adecuado para datos contables con sobredispersión, ya que permite que la varianza exceda a la media.