En esta práctica aprenderemos a ajustar el modelo GLM de Poisson en R a través de distintos ejemplos que consideran particularidades como el uso de offset, la sobredispersión, el exceso de ceros o la interacción entre predictores. En ellos se trabajará sobre la interpretación del modelo, su validez, inferencias sobre los parámetros y detección de puntos influyentes y/o atípicos. El objetivo es que el estudiante adquiera un dominio básico para abordar problemas similares en la práctica.

1. Ajuste del modelo, estimación, y residuos. Datos simulados

Generación de Datos

Simulamos datos de una variable de respuesta \(Y \sim Poisson(\lambda)\), donde:

\[ \log(\lambda) = \beta_0 + \beta_1 x \]

con \(\beta_0 = -1\) y \(\beta_1 = 0.3\), y donde la covariable \(x \sim U(0,10)\).

# Cargar librerías necesarias
library(ggplot2)
library(boot)
library(statmod)
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:boot':
## 
##     logit

library(glmtoolbox)

## 
## Attaching package: 'glmtoolbox'

## The following object is masked from 'package:boot':
## 
##     envelope

# Función para generar los datos
gen_dat <- function(n, b0, b1) {
  x <- runif(n=n, min=0, max=10)  # Covariable x ~ Uniforme(0,1)
  lambda <- exp(b0 + b1 * x)  # Relación log-lineal con Poisson
  y <- rpois(n=n, lambda=lambda)  # Respuesta Poisson
  data.frame(y=y, x=x)
}

# Generando datos para n = 150
n <- 150
datos <- gen_dat(n=n, b0=-1, b1=0.3)

# Visualizar los primeros datos
head(datos)

# Gráfico de dispersión de los datos simulados
ggplot(datos, aes(x = x, y = y)) +
  geom_point(alpha = 0.6, color = "blue") +  # Datos originales en azul
  labs(title = "Datos simulados: y ~ Poisson(lambda)",
       x = "Covariable x", 
       y = "Número de eventos (y)") +
  theme_minimal()

# Ajustamos un modelo GLM con familia Poisson
model_pois <- glm(y ~ x, family = poisson(link="log"), data = datos)

# Resumen del modelo
summary(model_pois)

## 
## Call:
## glm(formula = y ~ x, family = poisson(link = "log"), data = datos)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.95965    0.17304  -5.546 2.92e-08 ***
## x            0.29366    0.02162  13.581  < 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: 374.23  on 149  degrees of freedom
## Residual deviance: 130.34  on 148  degrees of freedom
## AIC: 477.98
## 
## Number of Fisher Scoring iterations: 5

# Extraemos los coeficientes del modelo
coef_ajustados <- coef(model_pois)
coef_exponenciales <- exp(coef_ajustados)  # Exponencial de los coeficientes

# Tabla de comparación
coef_comparacion <- data.frame(
  Parámetro = c("Intercept", "Pendiente"),
  Simulado = c(-1, 0.3),
  Estimado = coef_ajustados,
  `Exp(coef)` = coef_exponenciales
)
coef_comparacion

##             Parámetro Simulado   Estimado Exp.coef.
## (Intercept) Intercept     -1.0 -0.9596544 0.3830253
## x           Pendiente      0.3  0.2936562 1.3413226

# En un modelo de Poisson, los coeficientes NO representan diferencias absolutas, sino diferencias en el logaritmo de la tasa esperada de eventos. Para interpretarlos, aplicamos la función exponencial 

# El Intercept representa la tasa de incidencia esperada cuando x = 0.
#La pendiente (exp(beta1)) es la razón de tasas de incidencia a veces denominada IRR.
#Si beta1=0.29, entonces exp(0.3)=1.33, lo que significa que por cada unidad adicional en x, la tasa esperada de y aumenta en un 33%.

# Predicciones del modelo
datos$pred <- predict(model_pois, type = "response")

# Graficamos los valores observados vs. ajustados
ggplot(datos, aes(x = x, y = y)) +
  geom_point(alpha = 0.5, color = "blue") +
  geom_line(aes(y = pred), color = "red", size = 1) +
  labs(title = "Valores Observados vs. Ajustados",
       x = "Variable Predictora (x)",
       y = "Conteo de Eventos (y)",
       subtitle = "Los valores predichos (línea roja) siguen la tendencia de los datos") +
  theme_minimal()

# Coeficiente estimado positivo

Tarea. Genera un grid de tamaños mustrales (n) y evalua cómo el tamaño de muestra afecta a la estimación del de los coeficientes del modelo. Considera medidas como el valor medio de los coeficientes estimados o el error cuadrático medio, incluye gráficos de dispersión de estas medidas en función de los valores de n.

Residuos del modelo

# Calcular residuos
datos$res_pearson  <- boot::glm.diag(model_pois)$rp   # Residuos de Pearson
datos$res_deviance <- boot::glm.diag(model_pois)$rd   # Residuos de Deviance
datos$res_quantile <- statmod::qresid(model_pois)     # Residuos Cuantiles

# Ver primeras filas con residuos
head(datos)

##   y         x      pred res_pearson res_deviance res_quantile
## 1 1 0.2452403 0.4116269   0.9224386    0.7781787    0.9864178
## 2 7 8.1791573 4.2300154   1.3542939    1.2364234    1.1538829
## 3 0 2.8948295 0.8962196  -0.9523019   -1.3467582   -1.7057854
## 4 1 6.6612498 2.7086835  -1.0423785   -1.1982992   -1.2594186
## 5 8 9.3558075 5.9759090   0.8376106    0.7959680    0.6838174
## 6 1 1.9687596 0.6828255   0.3861674    0.3609048    0.2956270

# a) Dispersión de los residuos en función de x 
ggplot(datos, aes(x = x)) +
  geom_point(aes(y = res_pearson), color = "blue", alpha = 0.6) +
  geom_point(aes(y = res_deviance), color = "red", alpha = 0.6) +
  geom_point(aes(y = res_quantile), color = "green", alpha = 0.6) +
  labs(title = "Comparación de Residuos: Pearson (azul), Deviance (rojo) y Cuantiles (verde)",
       x = "Covariable x", y = "Residuos") +
  theme_minimal()

#Los residuos de Pearson  pueden ser más grandes en valores extremos.
#Los residuos de Deviance  suelen estar más ajustados en la región central.
#Los residuos Cuantiles  están más distribuidos en torno a 0 y deberían parecer normalizados.


# b) QQ plots para evaluar normalidad
par(mfrow = c(1, 3))  
qqPlot(datos$res_pearson, dist = "norm", mean = 0, sd = 1, main = "QQ-Plot: Residuos Pearson")

## [1] 128  76

qqPlot(datos$res_deviance, dist = "norm", mean = 0, sd = 1, main = "QQ-Plot: Residuos Deviance")

## [1] 128  92

qqPlot(datos$res_quantile, dist = "norm", mean = 0, sd = 1, main = "QQ-Plot: Residuos Cuantiles")

## [1] 128  92

#Residuos de Pearson y Deviance pueden mostrar más valores atípicos.
#Residuos Cuantiles deben seguir mejor la normalidad.


# c) Comparación de histogramas de residuos
par(mfrow = c(1, 3)) 
hist(datos$res_pearson, main = "Histograma: Residuos Pearson", col = "blue", breaks = 20)
hist(datos$res_deviance, main = "Histograma: Residuos Deviance", col = "red", breaks = 20)
hist(datos$res_quantile, main = "Histograma: Residuos Cuantiles", col = "green", breaks = 20)

#Residuos de Pearson pueden mostrar más variabilidad.
#Residuos de Deviance siguen la distribución de Poisson pero pueden ser sesgados.
#Residuos Cuantiles deberían parecer normales.

2. El uso del offset. Datos traffic_data

Analizar los accidentes de tráfico ocurridos en 30 ciudades a lo largo de un año, teniendo en cuenta el número de vehículos registrados en cada ciudad y el tipo de carretera predominante en la ciudad (urbana, autopista o secundaria).

Cargar los datos

getwd()

## [1] "C:/Users/Yolanda/OneDrive - Universidad de Valladolid/Docencia/MASTER25/Practicas/GLM"

load(file="C:/Users/Yolanda/OneDrive - Universidad de Valladolid/Docencia/MASTER25/Practicas/GLM/Datos/traffic_data.RData")

Visualizar los datos generados

head(traffic_data)

##       city vehicles  road_type lambda accidents
## 1 Ciudad_1   150913     Urbana    3.5         7
## 2 Ciudad_2   396270  Autopista    1.8         9
## 3 Ciudad_3   210399 Secundaria    2.5         5
## 4 Ciudad_4   442679  Autopista    1.8        10
## 5 Ciudad_5   470829     Urbana    3.5        17
## 6 Ciudad_6    32323 Secundaria    2.5         0

# Distribución de accidentes por tipo de carretera
ggplot(traffic_data, aes(x = road_type, y = accidents)) +
  geom_boxplot(fill = "skyblue", alpha = 0.6) +
  labs(title = "Distribución de accidentes por tipo de carretera",
       x = "Tipo de Carretera", y = "Número de Accidentes")

# Relación entre número de vehículos y número de accidentes 
ggplot(traffic_data, aes(x = vehicles, y = accidents, color = road_type)) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Número de accidentes  vs. Número de vehículos",
       x = "Número de Vehículos Registrados", y = "Número de Accidentes")

## `geom_smooth()` using formula = 'y ~ x'

Ajuste de un modelo de Poisson

# Modelo de Poisson sin offset
mod_no_offset <- glm(accidents ~ road_type, family = poisson, data = traffic_data)
# Nota: road_type se está tratando como un factor, dondo Autopista es la categoría de referencia. Orden alfabético

# Resumen del modelo
summary(mod_no_offset)

## 
## Call:
## glm(formula = accidents ~ road_type, family = poisson, data = traffic_data)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          1.66501    0.16440  10.128  < 2e-16 ***
## road_typeSecundaria -0.01635    0.21508  -0.076    0.939    
## road_typeUrbana      0.76048    0.18393   4.135 3.56e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 125.365  on 29  degrees of freedom
## Residual deviance:  91.166  on 27  degrees of freedom
## AIC: 199.26
## 
## Number of Fisher Scoring iterations: 5

#Este modelo no tiene en cuenta que algunas ciudades tienen más vehículos que otras, lo que puede sesgar la interpretación. No vamos a entrar en la interpretacaión de los coeficientes en este modelo

Incorporación offset

Un offset en un GLM es un término que se incluye en la ecuación del modelo cuya contribución al predictor lineal es fija y no se estima a partir de los datos. El uso de offsets es adecuado cuando la respuesta es un conteo relativo a un denominador conocido, como tiempo de exposición o número de unidades observadas. También puede resultar útil en situaciones donde una variable ya está en escala logarítmica y no debe ser estimada, en estos casos se usa como offset para mantener la relación adecuada.

El modelo estándar de Poisson es:

\[ \log(\lambda) = \beta_0 + \beta_1 x \]

Si queremos modelizar una tasa, necesitamos incluir un offset \( (t) \), donde \( t \) es el tiempo de exposición o el denominador de la tasa:

\[ \log(\lambda) = \beta_0 + \beta_1 x + \log(t). \]

Aquí, el término \( (t) \) no necesita un coeficiente estimado porque ya es un valor conocido, porque su efecto es fijo. Estamos modelizando la tasa \(\lambda/t\), y el término \(\log(t)\) simplemente ajusta la escala de la respuesta.

Ajuste de un modelo de Poisson con offset (ajustado por número de vehículos)

mod_offset <- glm(accidents ~ road_type + offset(log(vehicles)), 
                  family = poisson, 
                  data = traffic_data)

# Resumen del modelo
summary(mod_offset) # Cambio en los coeficientes del modelo

## 
## Call:
## glm(formula = accidents ~ road_type + offset(log(vehicles)), 
##     family = poisson, data = traffic_data)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -10.9614     0.1644 -66.676  < 2e-16 ***
## road_typeSecundaria   0.1815     0.2151   0.844    0.399    
## road_typeUrbana       0.7292     0.1839   3.964 7.36e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 53.315  on 29  degrees of freedom
## Residual deviance: 29.996  on 27  degrees of freedom
## AIC: 138.09
## 
## Number of Fisher Scoring iterations: 4

#El offset ajusta por la diferencia en el número de vehículos en cada ciudad

# log(vehicles) actúa como offset para modelizar la tasa de accidentes por vehículo en lugar del número absoluto de accidentes.

# En un modelo de Poisson, los coeficientes NO representan diferencias absolutas, sino diferencias en el logaritmo de la tasa esperada de eventos. Para interpretarlos, aplicamos la función exponencial 


# Para road_type=urbana (significativa):
#Interpretación en la escala logarítmica: El logaritmo de la tasa esperada de eventos es 0.73 unidades mayor que en la categoría de referencia (autopista)

#¿Interpretación en términos de tasa esperada?

Comparación de modelos

anova(mod_no_offset, mod_offset, test = "LRT") # ¿Qué ocurre? ¿Menor deviance?

## Analysis of Deviance Table
## 
## Model 1: accidents ~ road_type
## Model 2: accidents ~ road_type + offset(log(vehicles))
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        27     91.166                     
## 2        27     29.996  0    61.17

AIC(mod_no_offset, mod_offset)# ¿Menor AIC?

##               df      AIC
## mod_no_offset  3 199.2638
## mod_offset     3 138.0940

p_value <- pchisq ( deviance ( mod_offset ) , df.residual ( mod_offset ) , lower.tail = FALSE )
cat ( " Deviance : " , deviance ( mod_offset ) , " P-value : " , p_value  )

##  Deviance :  29.99593  P-value :  0.3143366

# El efecto del tipo de carretera  es estadísticamente significativo.
# El offset permite comparar las tasas de accidentes por vehículo, eliminando el sesgo de ciudades con mayor número de automóviles.

Predicciones del modelo con offset

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

traffic_data <- traffic_data %>%
  mutate(pred_accidents = predict(mod_offset, type = "response")) # crea nuevas variables o modifica variables existentes en un dataframe sin alterar las demás columnas.

# Gráfico de predicciones
ggplot(traffic_data, aes(x = road_type, y = pred_accidents)) +
  geom_boxplot(fill = "orange", alpha = 0.6) +
  labs(title = "Predicción del Modelo de Poisson con Offset",
       x = "Tipo de Carretera", y = "Tasa de Accidentes Ajustada")

3. Test de hipótesis y sobredispersión. Datos Crabs dataset de Agresti

Durante la temporada de desove, las hembras de una especie de cangrejo migran a la costa para realizar la puesta, con un macho enganchado en su cola. Cavan en la arena, donde depositan los huevos que son fertilizados externamente, tanto por el macho enganchado a la hembra como por otros machos que se reúnen alrededor de la pareja. Estos otros machos reciben el nombre de satélites. El archivo proporcionado contiene datos de un estudio realizado sobre esta especie en el golfo de Mexico.

Se desea modelizar el número de satélites de las cangrejas en función del peso, color, ancho de caparazón y el estado de la espina dorsal. En el modelo elegido, ¿se puede concluir que el peso del caparazón es significativamente mayor a 0.1? ¿Se cumplen las asunciones del modelo? ¿Son válidos estos resultados?

Carga de Datos, otra manera diferente

# Cargar datos desde URL
url <- "http://users.stat.ufl.edu/~aa/glm/data/Crabs.dat"
Crabs <- read.table(url, header=TRUE)
#weight: Peso del cangrejo (kg)
#width: Ancho del caparazón (cm)
#color: Clasificación del color (1: Claro, 4: Oscuro)
#spine: Condición de la espina dorsal (1: Buen estado, 3: Ambas desgastadas)
#y: number of satelites

# Convertir variables categóricas a factores
Crabs$color <- as.factor(Crabs$color)
Crabs$spine <- as.factor(Crabs$spine)

# Ver primeras filas
head(Crabs)

##   crab y weight width color spine
## 1    1 8   3.05  28.3     2     3
## 2    2 0   1.55  22.5     3     3
## 3    3 9   2.30  26.0     1     1
## 4    4 0   2.10  24.8     3     3
## 5    5 4   2.60  26.0     3     3
## 6    6 0   2.10  23.8     2     3

Ajuste del modelo GLM de Poisson con todas las covariables

fit_poisson <- glm(y ~ weight + color + width + spine, 
                   family = poisson(link="log"), data=Crabs)

# Resumen del modelo
summary(fit_poisson)

## 
## Call:
## glm(formula = y ~ weight + color + width + spine, family = poisson(link = "log"), 
##     data = Crabs)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -0.36180    0.96655  -0.374  0.70817   
## weight       0.49647    0.16626   2.986  0.00283 **
## color2      -0.26485    0.16811  -1.575  0.11515   
## color3      -0.51371    0.19536  -2.629  0.00855 **
## color4      -0.53086    0.22692  -2.339  0.01931 * 
## width        0.01675    0.04892   0.342  0.73207   
## spine2      -0.15037    0.21358  -0.704  0.48139   
## spine3       0.08728    0.11993   0.728  0.46674   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 632.79  on 172  degrees of freedom
## Residual deviance: 549.59  on 165  degrees of freedom
## AIC: 920.88
## 
## Number of Fisher Scoring iterations: 6

# width y spine no resultan significativas
# ¿Cuántas observaciones hay? ¿Qué relación hay con el número de gl de residual y null deviance?
# ¿Qué indica 'Number of Fisher Scoring Iterations'?

Ajuste del modelo GLM de Poisson con covariables significativas

fit_poisson2 <- glm(y ~ weight + color , 
                   family = poisson(link="log"), data=Crabs)

# Resumen del modelo
summary(fit_poisson2)

## 
## Call:
## glm(formula = y ~ weight + color, family = poisson(link = "log"), 
##     data = Crabs)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.04978    0.23315  -0.214   0.8309    
## weight       0.54618    0.06811   8.019 1.07e-15 ***
## color2      -0.20511    0.15371  -1.334   0.1821    
## color3      -0.44980    0.17574  -2.560   0.0105 *  
## color4      -0.45205    0.20844  -2.169   0.0301 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 632.79  on 172  degrees of freedom
## Residual deviance: 551.80  on 168  degrees of freedom
## AIC: 917.1
## 
## Number of Fisher Scoring iterations: 6

Interpretacino de los coeficientes

# En un modelo de Poisson, los coeficientes NO representan diferencias absolutas, sino diferencias en el logaritmo de la tasa esperada de eventos. Para interpretarlos, aplicamos la función exponencial (tasa esperada)

exp(coef(fit_poisson2))

## (Intercept)      weight      color2      color3      color4 
##   0.9514381   1.7266442   0.8145557   0.6377580   0.6363225

exp(confint(fit_poisson2)) # Cálculo de intervalos de confianza ¿Cómo? ¿Que información contiene cada columna de la tabla de coeficientes estimados?

##                 2.5 %    97.5 %
## (Intercept) 0.6009718 1.5001076
## weight      1.5076222 1.9689147
## color2      0.6089529 1.1138967
## color3      0.4544755 0.9066422
## color4      0.4218499 0.9575813

# Por cada unidad adicional de peso, la tasa esperada de eventos aumenta en un (exp(0.55)) 73%, manteniendo constante el color.

# Para color 3:
#Interpretación en la escala logarítmica:  el logaritmo de la tasa esperada de eventos es 0.45 unidades menor que para los individuos con Color1 (categoría de referencia), manteniendo constante el peso.
#Interpretación en términos de tasa esperada: Aplicamos la exponencial exp(-0.45)=0.64. Esto significa que los individuos con Color3 tienen una tasa esperada de eventos un 36% menor (1-0.64) en comparación con los de Color1, manteniendo constante el peso.

Gráfico de predicciones

Crabs$predicted <- predict(fit_poisson2, type="response")

# Gráfico con jitter para visualizar mejor las coincidencias
ggplot(Crabs, aes(x = weight, y = predicted, color = color)) +
  geom_point(aes(y = y)) +
  geom_line(size = 1) +
  labs(title = "Predicción del Número de Satélites",
       x = "Peso",
       y = "Número Esperado de Satelite") +
  theme_minimal()

Test de Wald para coeficientes

# Calculando el estadistico y su valor-P para contrastar H0: beta_1=0.1; H1: beta_1>0.1
z <- (0.546 - 0.1) / 0.068
z

## [1] 6.558824

pnorm(q=z, lower.tail=FALSE) # valor-P

## [1] 2.711697e-11

# El p-valor es menor grande que un nivel de significanción usual, se rechaza H0

Bondad de ajuste

# Test de bondad de ajuste para el modelo de Poisson
pchisq(fit_poisson2$deviance, df=fit_poisson2$df.residual, lower.tail=FALSE)

## [1] 2.08043e-42

# Posibles causas: falta de linealidad en el logaritmo de la respuesta, exceso de variabilidad en la misma, o quizás ausencia en el modelo de variables importantes que debieran tenerse en cuenta (y que el investigador ignora)

Sobredispersión?

# 1. Comparación entre media y varianza
mean_sat <- mean(Crabs$y)
var_sat <- var(Crabs$y)
cat("Media de premios:", mean_sat, " - Varianza de premios:", var_sat, "\n")

## Media de premios: 2.919075  - Varianza de premios: 9.912018

#Consecuencias de la Sobredispersión:
# Estimaciones Sesgadas: En presencia de sobredispersión, las estimaciones de los parámetros del modelo de Poisson pueden estar sesgadas.

# Errores Estándar Subestimados: Los errores estándar de las estimaciones pueden subestimarse, llevando a inferencias estadísticas incorrectas, como intervalos de confianza demasiado estrechos o pruebas de hipótesis que rechazan incorrectamente la hipótesis nula.

# 2. Relación entre deviance y grados de libertad
deviance_ratio <- deviance(fit_poisson2) / df.residual(fit_poisson2)
cat("Razón Deviance/GL:", deviance_ratio, "\n")

## Razón Deviance/GL: 3.284553

if (deviance_ratio > 1.5) {
  cat("Indicación de sobre-dispersión en el modelo.\n")
} else {
  cat("No hay evidencia fuerte de sobre-dispersión.\n")
}

## Indicación de sobre-dispersión en el modelo.

# 3. Test de sobre-dispersión de Cameron & Trivedi. Para GLM de Poisson 
# H1: Var > media
# Estadístico basado en los residuos de Pearson, bajo H0 chi1 con n-p
pearson_chisq <- sum(residuals(fit_poisson2, type = "pearson")^2)
dispersion_ratio <- pearson_chisq / df.residual(fit_poisson2)
p_value_dispersion <- pchisq(pearson_chisq, df.residual(fit_poisson2), lower.tail=FALSE)

cat("Pearson Chi2:", pearson_chisq, " - p-value:", p_value_dispersion, "\n")

## Pearson Chi2: 534.8636  - p-value: 7.564098e-40

if (p_value_dispersion < 0.05) {
  cat("El test de sobre-dispersión indica que el modelo de Poisson puede no ser adecuado.\n")
}

## El test de sobre-dispersión indica que el modelo de Poisson puede no ser adecuado.

# 4. Paquete espefífico test de sobredispersión
# Test de sobredispersión
#library(AER)
#dispersiontest(modPois2)

# Cuando se detecta sobredispersión en un modelo de Poisson, se puede utilizar como alternativa la regresión binomial negativa.

¿Son válidas las inferencias e interpretación del modelo que hemos hecho?

4. Selección modelo, evaluación y diagnóstico. Datos awards

En una escuela secundaria, se ha recopilado información sobre el número de premios académicos que han obtenido los estudiantes en un año. Se quiere analizar qué factores influyen en la cantidad de premios ganados y si ciertas características de los estudiantes pueden predecir su desempeño en la obtención de premios.

Cargar los datos

# Cargar datos
awards <- read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")
awards <- as.data.frame(awards)

# Convertir la variable de programa a factor
awards <- within(awards, {
  prog <- factor(prog, levels = 1:3, labels = c("General", "Academic", "Vocational"))
})

Resumen de los datos

summary(awards)

##        id           num_awards           prog          math      
##  Min.   :  1.00   Min.   :0.00   General   : 45   Min.   :33.00  
##  1st Qu.: 50.75   1st Qu.:0.00   Academic  :105   1st Qu.:45.00  
##  Median :100.50   Median :0.00   Vocational: 50   Median :52.00  
##  Mean   :100.50   Mean   :0.63                    Mean   :52.65  
##  3rd Qu.:150.25   3rd Qu.:1.00                    3rd Qu.:59.00  
##  Max.   :200.00   Max.   :6.00                    Max.   :75.00

ggplot(awards, aes(x = math, y = num_awards, color = prog)) +
  geom_point(alpha = 0.7, size = 3) +
  labs(title = "Número de Premios en función del Puntaje de Matemáticas",
       x = "Puntaje en Matemáticas",
       y = "Número de Premios",
       color = "Programa") +
  theme_minimal() +
  theme(legend.position = "top")

Asunciones GLM Poisson

# Media y desviación estándar del número de premios
mean_awards <- mean(awards$num_awards)
sd_awards <- sd(awards$num_awards)
cat("Media de premios:", mean_awards, " - Desviación estándar:", sd_awards, "\n")

## Media de premios: 0.63  - Desviación estándar: 1.052921

# Media y desviación estándar por tipo de programa
with(awards, tapply(num_awards, prog, function(x) {
  sprintf("M (SD) = %1.2f (%1.2f)", mean(x), sd(x))
}))

##                General               Academic             Vocational 
## "M (SD) = 0.20 (0.40)" "M (SD) = 1.00 (1.28)" "M (SD) = 0.24 (0.52)"

Ajuste de modelos anidados

# Ajuste de modelos GLM de Poisson
fit_glm_0 <- glm(num_awards ~ 1, data=awards, family=poisson)
fit_glm_1 <- glm(num_awards ~ math, data=awards, family=poisson)
fit_glm_2 <- glm(num_awards ~ math + prog, data=awards, family=poisson)

# Resumen de los modelos
summary(fit_glm_0)

## 
## Call:
## glm(formula = num_awards ~ 1, family = poisson, data = awards)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.46204    0.08909  -5.186 2.14e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 287.67  on 199  degrees of freedom
## Residual deviance: 287.67  on 199  degrees of freedom
## AIC: 465.73
## 
## Number of Fisher Scoring iterations: 6

summary(fit_glm_1)

## 
## Call:
## glm(formula = num_awards ~ math, family = poisson, data = awards)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.333532   0.591261  -9.021   <2e-16 ***
## math         0.086166   0.009679   8.902   <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: 287.67  on 199  degrees of freedom
## Residual deviance: 204.02  on 198  degrees of freedom
## AIC: 384.08
## 
## Number of Fisher Scoring iterations: 6

summary(fit_glm_2)

## 
## Call:
## glm(formula = num_awards ~ math + prog, family = poisson, data = awards)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -5.24712    0.65845  -7.969 1.60e-15 ***
## math            0.07015    0.01060   6.619 3.63e-11 ***
## progAcademic    1.08386    0.35825   3.025  0.00248 ** 
## progVocational  0.36981    0.44107   0.838  0.40179    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 287.67  on 199  degrees of freedom
## Residual deviance: 189.45  on 196  degrees of freedom
## AIC: 373.5
## 
## Number of Fisher Scoring iterations: 6

Tarea. Interpretación del modelo fit_glm_2.

Sobredispersión?

# 1. Comparación entre media y varianza
mean_awards <- mean(awards$num_awards)
var_awards <- var(awards$num_awards)
cat("Media de premios:", mean_awards, " - Varianza de premios:", var_awards, "\n")

## Media de premios: 0.63  - Varianza de premios: 1.108643

# 2. Relación entre deviance y grados de libertad
deviance_ratio <- deviance(fit_glm_2) / df.residual(fit_glm_2)
cat("Razón Deviance/GL:", deviance_ratio, "\n")

## Razón Deviance/GL: 0.9665797

if (deviance_ratio > 1.5) {
  cat("Indicación de sobre-dispersión en el modelo.\n")
} else {
  cat("No hay evidencia fuerte de sobre-dispersión.\n")
}

## No hay evidencia fuerte de sobre-dispersión.

# 3. Test de sobre-dispersión de Cameron & Trivedi
pearson_chisq <- sum(residuals(fit_glm_2, type = "pearson")^2)
dispersion_ratio <- pearson_chisq / df.residual(fit_glm_2)
p_value_dispersion <- pchisq(pearson_chisq, df.residual(fit_glm_2), lower.tail=FALSE)

cat("Pearson ChiÃÂ²:", pearson_chisq, " - p-value:", p_value_dispersion, "\n")

## Pearson ChiÃÂ²: 212.1437  - p-value: 0.203986

if (p_value_dispersion < 0.05) {
  cat("El test de sobre-dispersión indica que el modelo de Poisson puede no ser adecuado.\n")
}

# 4. Librería
# Test de sobredispersión
#library(AER)
#dispersiontest(modPois2)

Comparacion y seleccion de modelos. Test LRT y Deviance

# Comparación de modelos con LRT
anova(fit_glm_1, fit_glm_2, test="LRT")# Tamaño muestral moderado-grande convergencia asintótica a chi2

## Analysis of Deviance Table
## 
## Model 1: num_awards ~ math
## Model 2: num_awards ~ math + prog
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       198     204.02                          
## 2       196     189.45  2   14.572 0.0006852 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Bondad de ajuste

# Prueba de bondad de ajuste usando deviance
p_value <- pchisq(deviance(fit_glm_2), df.residual(fit_glm_2), lower.tail=FALSE)
cat("Deviance:", deviance(fit_glm_2), " - p-value:", p_value, "\n")

## Deviance: 189.4496  - p-value: 0.6182274

# Cálculo de pseudo R2 de McFadden
logLik_full <- logLik(fit_glm_2)
logLik_null <- logLik(fit_glm_0)
pseudo_R2 <- 1 - (logLik_full / logLik_null)
cat("Pseudo R2 de McFadden:", pseudo_R2, "\n")

## Pseudo R2 de McFadden: 0.2118112

# Otras maneras de calcular R2
#library ( DescTools )
#PseudoR2(fit_glm_2 , which = "all" )

Análisis de residuos

par(mfrow=c(2,2)) # Configurar panel de gráficos
plot(fit_glm_2) # Gráficos de diagnóstico estándar

# Extraer residuos deviance
residuos_dev <- residuals(fit_glm_2, type="deviance")

# Histogramas y gráficos de dispersión
hist(residuos_dev, main="Histograma de Residuos Deviance", col="lightblue", breaks=20)
qqnorm(residuos_dev)
qqline(residuos_dev, col="red", lwd=2)

# Identificación de valores atípicos
outliers <- which(abs(residuos_dev) > 2) # Considerar valores fuera de ÃÂ±2
cat("Observaciones con residuos altos:", outliers, "\n")

## Observaciones con residuos altos: 54 157 164 181 199

Diagnostico de puntos influyentes y outliers

# Cálculo de medidas de influencia
cooks_d <- cooks.distance(fit_glm_2)
hat_values <- hatvalues(fit_glm_2)

# Identificación de puntos influyentes
threshold_cooks <- 4 / (nrow(awards) - length(coef(fit_glm_2)))
outliers <- which(cooks_d > threshold_cooks)
outliers

##  27  28  36  54  58  87 154 157 181 190 194 199 
##  27  28  36  54  58  87 154 157 181 190 194 199

cat("Puntos influyentes identificados:", outliers, "\n")

## Puntos influyentes identificados: 27 28 36 54 58 87 154 157 181 190 194 199

# Visualización de CookÃ¢ÂÂs Distance
plot(cooks_d, type="h", col="red", lwd=2, ylab="Cook's Distance", xlab="Observación", 
     main="Distancia de Cook")
abline(h=threshold_cooks, col="blue", lty=2)

# Eliminación de outliers y ajuste de nuevo modelo
awards_clean <- awards[-outliers, ]
fit_glm_2_clean <- glm(num_awards ~ math + prog, data=awards_clean, family=poisson)
summary(fit_glm_2_clean)

## 
## Call:
## glm(formula = num_awards ~ math + prog, family = poisson, data = awards_clean)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -6.26324    0.80904  -7.742 9.82e-15 ***
## math            0.07913    0.01239   6.388 1.68e-10 ***
## progAcademic    1.45359    0.46829   3.104  0.00191 ** 
## progVocational  0.66336    0.55782   1.189  0.23437    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 235.17  on 187  degrees of freedom
## Residual deviance: 143.19  on 184  degrees of freedom
## AIC: 299.02
## 
## Number of Fisher Scoring iterations: 6

Tarea. Compara los resultados del modelo fit_glm_2_clean con fit_glm_2.

Predicciones con el modelo sin outliers

awards_clean$predicted <- predict(fit_glm_2_clean, type="response")

# Gráfico con jitter para visualizar mejor las coincidencias
ggplot(awards_clean, aes(x = math, y = predicted, color = prog)) +
  geom_point(aes(y = num_awards), alpha = 0.5, position = position_jitter(h = 0.2)) +
  geom_line(size = 1) +
  labs(title = "Predicción del Número de Premios sin Outliers",
       x = "Puntaje en Matemáticas",
       y = "Número Esperado de Premios") +
  theme_minimal()

5. GLM Poison con ceros inflados. Datos fish

Se recopilaron datos sobre 250 grupos de visitantes que acudieron a un parque estatal donde se permite la pesca. A los visitantes se les realizaron algunas preguntas y se obtuvo información sobre las siguientes variables:

persons: Número de personas en el grupo.
child: Número de niños en el grupo.
camper: Indicador de si el grupo llevó una caravana al parque o no.
count: Número de peces capturados por el grupo.

El objetivo del estudio es predecir el número de peces capturados por un grupo en función de su composición. Dado que la variable de interés es una variable aleatoria discreta que solo toma valores enteros no negativos, se opta por ajustar un modelo de regresión para datos de conteo.

Al analizar los datos, observamos una cantidad excesiva de ceros, lo cual puede deberse a dos razones principales:

Algunos grupos no pescaron en absoluto, pero como no tenemos una variable que registre si la actividad de pesca ocurrió o no, estas observaciones se registraron como ceros.
Otros grupos intentaron pescar pero no lograron capturar ningún pez, lo que también resultó en un conteo de cero.

Este fenómeno genera una varianza mucho mayor en los datos en comparación con su media, lo que viola una de las suposiciones clave del modelo de regresión de Poisson.

Debido a esta sobre-dispersión y el exceso de ceros en los datos, se considera la aplicación de modelos alternativos, tales como:

Modelo binomial negativo: Permite modelizar datos con mayor varianza que la esperada en Poisson.
Modelos inflados en ceros (Zero-Inflated Poisson, ZIP): Diferencian entre ceros estructurales (grupos que no pescaron) y ceros aleatorios (grupos que pescaron pero no capturaron peces).
Modelos hurdle: Separan el proceso de decisión de pescar del proceso de captura de peces.

data <-
  within(read.csv("C:/Users/Yolanda/OneDrive - Universidad de Valladolid/Docencia/MASTER25/Practicas/GLM/Datos/fish.csv"),
         camper <- factor(camper))
attach(data)


######  Descriptive Analysis    ######

n <- length(count)
xi <- sort(unique(count))
fi <- table(count)

am <- mean(count)
variance <-var(count)

print(c(am, variance))

## [1]   3.2960 135.3739

plot(
  xi,
  fi,
  type = "h",
  lwd = 2,
  xlab = "No. of Fishes Caught",
  ylab = "Frequency",
  main = "Column Diagram of Fishing data"
)

6. Interacción, VIF y colinealidad. Datos gotelli

Gotelli & Ellison (2002) analizaron los determinantes biogeográficos de la riqueza de hormigas (Srich) a escala regional (gotelli.csv). Para esto se describieron el tipo de hábitat (Habitat), la latitud (Latitude) y la altitud (Elevation). Proponer un modelo adecuado para el análisis de la riqueza de hormigas a partir de las variables disponibles.

Carga de datos y exploración inicial

# Cargar paquetes
if (!require(car)){install.packages("car", dependencies=TRUE)}
if (!require(MuMIn)){install.packages("MuMIn", dependencies=TRUE)}

## Loading required package: MuMIn

## 
## Attaching package: 'MuMIn'

## The following object is masked from 'package:glmtoolbox':
## 
##     QIC

library(car)
library(MuMIn)

# Cargar datos
gotelli <- read.csv("C:/Users/Yolanda/OneDrive - Universidad de Valladolid/Docencia/MASTER25/Practicas/GLM/Datos/gotelli.csv")  
head(gotelli)

##   Site Srich Habitat Latitude Elevation
## 1  TPB     6  Forest    41.97       389
## 2  HBC    16  Forest    42.00         8
## 3  CKB    18  Forest    42.03       152
## 4  SKP    17  Forest    42.05         1
## 5   CB     9  Forest    42.05       210
## 6   RP    15  Forest    42.17        78

# Convertir la variable categórica Habitat a numérica antes de scatterplotMatrix
gotelli$HabitatNum <- as.numeric(as.factor(gotelli$Habitat))

# Visualización inicial
# Crear la matriz de dispersión con la nueva variable numérica
scatterplotMatrix(~ Srich + HabitatNum + Latitude + Elevation, data = gotelli)

pairs(gotelli[,c(2,4,5,6)])

Asunciones del modelo

# Ajuste de distribución a los datos
sim.pois <- dpois(x = 0:max(gotelli$Srich), lambda = mean(gotelli$Srich))
plot(table(gotelli$Srich), xlab = "Número de especies", ylab = "Frecuencia")

hist(gotelli$Srich, xlab = "Número de especies", ylab = "Frecuencia relativa", main = "", freq = FALSE)
lines(x = 0:max(gotelli$Srich), y = sim.pois, col = "blue", lwd = 2, type = "b")

var(gotelli$Srich)/mean(gotelli$Srich)

## [1] 2.566343

Ajuste del modelo de Poison inicial con todas interacciones

# Ajustar el modelo GLM de Poisson
gotelli.glm <- glm(Srich ~ Habitat * Latitude * Elevation, family = poisson, data = gotelli)
summary(gotelli.glm)

## 
## Call:
## glm(formula = Srich ~ Habitat * Latitude * Elevation, family = poisson, 
##     data = gotelli)
## 
## Coefficients:
##                                    Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                      16.1115492  8.2821982   1.945   0.0517 .
## HabitatForest                    -2.1996735 10.0663651  -0.219   0.8270  
## Latitude                         -0.3355369  0.1950119  -1.721   0.0853 .
## Elevation                        -0.0209554  0.0307527  -0.681   0.4956  
## HabitatForest:Latitude            0.0690150  0.2369897   0.291   0.7709  
## HabitatForest:Elevation           0.0137995  0.0381499   0.362   0.7176  
## Latitude:Elevation                0.0004718  0.0007208   0.655   0.5127  
## HabitatForest:Latitude:Elevation -0.0003348  0.0008941  -0.375   0.7080  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 102.763  on 43  degrees of freedom
## Residual deviance:  39.772  on 36  degrees of freedom
## AIC: 216.13
## 
## Number of Fisher Scoring iterations: 4

# Evaluación de colinealidad mediante VIF (Factor de Inflación de la varianza): métrica utilizada para detectar la multicolinealidad en modelos de regresión. La multicolinealidad ocurre cuando dos o más variables explicativas están altamente correlacionadas, lo que puede distorsionar la estimación de los coeficientes y afectar la interpretabilidad del modelo.
vif_values <- vif(gotelli.glm)

## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif

print(vif_values)  # Hay problemas de colinealidad si los valores son altos (>10) indican multicolinealidad severa. Se requieren  < 5

##                    Habitat                   Latitude 
##                 7088.06105                   10.55347 
##                  Elevation           Habitat:Latitude 
##                 7124.53754                 7184.92459 
##          Habitat:Elevation         Latitude:Elevation 
##                10812.59771                 7247.74145 
## Habitat:Latitude:Elevation 
##                10971.71840

#  Centrar las variables (restar la media) puede ayudar a reducir la colinealidad en algunos casos, especialmente cuando hay interacciones entre variables.

Reducir colinealidad. Centrar la latitud y la elevación

gotelli$cLatitude <- scale(gotelli$Latitude, scale = TRUE)
gotelli$cElevation <- scale(gotelli$Elevation, scale = TRUE)

Reajuste del modelo con variables centradas

# Ajustar el modelo con variables centradas
gotelli.glm <- glm(Srich ~ Habitat * cLatitude * cElevation, family = poisson, data = gotelli,na.action = na.fail)
summary(gotelli.glm)

## 
## Call:
## glm(formula = Srich ~ Habitat * cLatitude * cElevation, family = poisson, 
##     data = gotelli, na.action = na.fail)
## 
## Coefficients:
##                                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                         1.523727   0.104428  14.591  < 2e-16 ***
## HabitatForest                       0.628476   0.129210   4.864 1.15e-06 ***
## cLatitude                          -0.241525   0.113340  -2.131   0.0331 *  
## cElevation                         -0.105970   0.110845  -0.956   0.3391    
## HabitatForest:cLatitude            -0.009535   0.140664  -0.068   0.9460    
## HabitatForest:cElevation           -0.097558   0.137497  -0.710   0.4780    
## cLatitude:cElevation                0.081362   0.124289   0.655   0.5127    
## HabitatForest:cLatitude:cElevation -0.057740   0.154172  -0.375   0.7080    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 102.763  on 43  degrees of freedom
## Residual deviance:  39.772  on 36  degrees of freedom
## AIC: 216.13
## 
## Number of Fisher Scoring iterations: 4

# Evaluar nuevamente la colinealidad
vif_values_new <- vif(gotelli.glm)

## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif

print(vif_values_new)

##                      Habitat                    cLatitude 
##                     1.167807                     3.113812 
##                   cElevation            Habitat:cLatitude 
##                     3.563564                     3.220434 
##           Habitat:cElevation         cLatitude:cElevation 
##                     3.609016                     3.477485 
## Habitat:cLatitude:cElevation 
##                     3.644151

# La colinealidad ha disminuido considerablemente y el modelo es más estable.

# Ayuda a la interpretación en el caso de colinealidad:
# HabitatForest:Mide el efecto de estar en un hábitat forestal en comparación con un hábitat de tipo Bog (referencia), cuando Latitud y Elevación son cero. El coeficiente positvo indica que, en promedio, la riqueza de especies en un bosque es mayor que en un bog, y este efecto no es significativo.
# Latitude: Mide el efecto de la latitud en la riqueza de especies, cuando Elevación = 0 y el hábitat es Bog. Un coeficiente negativo sugiere que a mayor latitud, menor riqueza de especies, aunque el efecto es marginalmente significativo.
# HabitatForest:Elevation Representa la modificación del efecto de la elevación en la riqueza de especies cuando el hábitat es un bosque. No significativo
# Latitude:Elevation: Representa la interacción entre latitud y elevación, es decir, si el efecto de la latitud cambia según la elevación. No significativo
# HabitatForest:Latitude:Elevation Representa el efecto de la latitud y la elevación en la riqueza de especies en función del hábitat. No significativo

Evaluación de la sobre-dispersión

dispersion_ratio <- gotelli.glm$deviance / gotelli.glm$df.residual
cat("Índice de sobre-dispersión:", dispersion_ratio, "\n")

## Índice de sobre-dispersión: 1.104765

# Si el índice es significativamente mayor que 1, el modelo de Poisson no es adecuado.

Comparación de modelos con AIC

summary(model.avg(dredge(gotelli.glm), fit = TRUE, subset = TRUE)) #  Explora todos los modelos posibles con combinaciones de predictores a partir del modelo gotelli.glm.

## Fixed term is "(Intercept)"
## Fixed term is "(Intercept)"
## Fixed term is "(Intercept)"

## 
## Call:
## model.avg(object = get.models(object = dredge(gotelli.glm), subset = TRUE))
## 
## Component model call: 
## glm(formula = Srich ~ <19 unique rhs>, family = poisson, data = 
##      gotelli, na.action = na.fail)
## 
## Component models: 
##         df  logLik   AICc delta weight
## 123      4 -100.52 210.07  0.00   0.43
## 1235     5 -100.31 212.21  2.14   0.15
## 1234     5 -100.34 212.26  2.19   0.14
## 1236     5 -100.50 212.58  2.51   0.12
## 12345    6 -100.13 214.54  4.47   0.05
## 12356    6 -100.31 214.89  4.82   0.04
## 12346    6 -100.32 214.91  4.84   0.04
## 23       3 -105.32 217.24  7.17   0.01
## 123456   7 -100.13 217.38  7.31   0.01
## 236      4 -105.30 219.62  9.55   0.00
## 1234567  8 -100.06 220.24 10.17   0.00
## 13       3 -108.50 223.60 13.53   0.00
## 135      4 -108.28 225.59 15.52   0.00
## 12       3 -115.37 237.33 27.26   0.00
## 3        2 -116.72 237.72 27.65   0.00
## 124      4 -115.18 239.40 29.33   0.00
## 2        2 -120.16 244.62 34.55   0.00
## 1        2 -123.34 250.98 40.91   0.00
## (Null)   1 -131.56 265.21 55.14   0.00
## 
## Term codes: 
##                   cElevation                    cLatitude 
##                            1                            2 
##                      Habitat         cElevation:cLatitude 
##                            3                            4 
##           cElevation:Habitat            cLatitude:Habitat 
##                            5                            6 
## cElevation:cLatitude:Habitat 
##                            7 
## 
## Model-averaged coefficients:  
## (full average) 
##                                      Estimate Std. Error Adjusted SE z value
## (Intercept)                         1.5292792  0.0995464   0.1026987  14.891
## cElevation                         -0.1647961  0.0792477   0.0813992   2.025
## cLatitude                          -0.2483057  0.0785800   0.0810422   3.064
## HabitatForest                       0.6300487  0.1214245   0.1252695   5.030
## cElevation:HabitatForest           -0.0196356  0.0704268   0.0721761   0.272
## cElevation:cLatitude                0.0107758  0.0412472   0.0423015   0.255
## cLatitude:HabitatForest            -0.0048191  0.0639284   0.0659522   0.073
## cElevation:cLatitude:HabitatForest -0.0001542  0.0085061   0.0087661   0.018
##                                    Pr(>|z|)    
## (Intercept)                         < 2e-16 ***
## cElevation                          0.04291 *  
## cLatitude                           0.00218 ** 
## HabitatForest                         5e-07 ***
## cElevation:HabitatForest            0.78558    
## cElevation:cLatitude                0.79893    
## cLatitude:HabitatForest             0.94175    
## cElevation:cLatitude:HabitatForest  0.98597    
##  
## (conditional average) 
##                                    Estimate Std. Error Adjusted SE z value
## (Intercept)                         1.52928    0.09955     0.10270  14.891
## cElevation                         -0.16741    0.07709     0.07933   2.110
## cLatitude                          -0.24847    0.07834     0.08081   3.075
## HabitatForest                       0.63005    0.12142     0.12527   5.030
## cElevation:HabitatForest           -0.07942    0.12375     0.12776   0.622
## cElevation:cLatitude                0.04444    0.07430     0.07671   0.579
## cLatitude:HabitatForest            -0.02218    0.13575     0.14013   0.158
## cElevation:cLatitude:HabitatForest -0.05774    0.15417     0.15953   0.362
##                                    Pr(>|z|)    
## (Intercept)                         < 2e-16 ***
## cElevation                          0.03483 *  
## cLatitude                           0.00211 ** 
## HabitatForest                         5e-07 ***
## cElevation:HabitatForest            0.53418    
## cElevation:cLatitude                0.56233    
## cLatitude:HabitatForest             0.87422    
## cElevation:cLatitude:HabitatForest  0.71740    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# En el Full Model Average, los coeficientes tienden a ser más pequeños porque las variables que no aparecen en ciertos modelos reciben coeficientes de 0.

# En el Conditional Model Average, los coeficientes solo se calculan para los modelos donde la variable está presente, por lo que suelen ser más grandes.

Tarea. A la vista de los resultados ajustar el modelo más adecuado, interpretar los resultados, evaluar la adecuación del modelo ási como la presencia de valores outliers o influyentes.

GLM de Poisson en R

Yolanda Larriba

2025-02-24

1. Ajuste del modelo, estimación, y residuos. Datos simulados

Generación de Datos

Residuos del modelo

2. El uso del offset. Datos traffic_data

Cargar los datos

Visualizar los datos generados

Ajuste de un modelo de Poisson

Incorporación offset

Ajuste de un modelo de Poisson con offset (ajustado por número de vehículos)

Comparación de modelos

Predicciones del modelo con offset

3. Test de hipótesis y sobredispersión. Datos Crabs dataset de Agresti

Carga de Datos, otra manera diferente

Ajuste del modelo GLM de Poisson con todas las covariables

Ajuste del modelo GLM de Poisson con covariables significativas

Interpretacino de los coeficientes

Gráfico de predicciones

Test de Wald para coeficientes

Bondad de ajuste

Sobredispersión?

4. Selección modelo, evaluación y diagnóstico. Datos awards

Cargar los datos

Resumen de los datos

Asunciones GLM Poisson

Ajuste de modelos anidados

Sobredispersión?

Comparacion y seleccion de modelos. Test LRT y Deviance

Bondad de ajuste

Análisis de residuos

Diagnostico de puntos influyentes y outliers

Predicciones con el modelo sin outliers

5. GLM Poison con ceros inflados. Datos fish

6. Interacción, VIF y colinealidad. Datos gotelli

Carga de datos y exploración inicial

Asunciones del modelo

Ajuste del modelo de Poison inicial con todas interacciones

Reducir colinealidad. Centrar la latitud y la elevación

Reajuste del modelo con variables centradas

Evaluación de la sobre-dispersión

Comparación de modelos con AIC