Modelos Lineales Generalizados: Análisis de datos provenientes de un ensayo en un DCL y dosis-respuesta
Olinda, R. A.
Noviembre, 2020
Ricardo Alves de Olinda
ricardo.estat@yahoo.com.br
http://lattes.cnpq.br/7767223263366578
Universidad del Estado de Paraíba
http://departamentos.uepb.edu.br/estatistica/corpo-docente/
Universidad de San Carlos de Guatemala
Centro de Telemática (CETE)
Facultad de Agronomía USAC
Use R!
R Markdown
Sobre Programación básica en R, puede ser consultada Aqui
Ejemplo 1
Los datos de este ejemplo se refieren a una serie de tres experimentos en diseño cuadrado latino 4X4. Los tratamientos son atrayentes sexuales para machos. Se contabilizó el número de insectos capturados cada semana, durante 8 semanas (Marino Barrientos, 2018).
Comienza cargando las librerías y preparando los datos a utilizar.
library(readxl)## Warning: package 'readxl' was built under R version 3.6.3library(MASS)## Warning: package 'MASS' was built under R version 3.6.3library(lattice)## Warning: package 'lattice' was built under R version 3.6.3library(latticeExtra)
library(multcomp)## Warning: package 'multcomp' was built under R version 3.6.3## Loading required package: mvtnorm## Warning: package 'mvtnorm' was built under R version 3.6.3## Loading required package: survival## Warning: package 'survival' was built under R version 3.6.3## Loading required package: TH.data##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserlibrary(hnp)## Warning: package 'hnp' was built under R version 3.6.3library(agricolae)## Warning: package 'agricolae' was built under R version 3.6.3library(htmltools)## Warning: package 'htmltools' was built under R version 3.6.3require(car)## Loading required package: car## Warning: package 'car' was built under R version 3.6.3## Loading required package: carData## Warning: package 'carData' was built under R version 3.6.3
La fotografía muestra una ventana en el pasillo de la Universidad de Cambridge, Caius College. Erigido en homenaje a Ronald Aylmer Fisher (cuadrado latino 7 x 7) y John Venn (diagrama de Venn)
Para realizar los análisis exploratorios a seguir, necesitaremos crear un objeto del tipo dataframe, como sigue.
#Importacion del datos para el R
da <- read_excel("Datos_conteo.xlsx",
sheet = "Data_R1")
head(da)## # A tibble: 6 x 7
## Experimento Semana Fila Columna Tratamiento Silvestres Laboratorio
## <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 1 1 1 2 A 0 0
## 2 1 2 1 2 A 0 0
## 3 1 3 1 2 A 1 0
## 4 1 4 1 2 A 1 0
## 5 1 5 1 2 A 0 0
## 6 1 6 1 2 A 0 0dim(da)## [1] 384 7# Nombre de columnas con las tres primeras letras
names(da) <- substr(tolower(names(da)), 0, 3)
names(da)[ncol(da)] <- "lab"
str(da)## tibble [384 x 7] (S3: tbl_df/tbl/data.frame)
## $ exp: num [1:384] 1 1 1 1 1 1 1 1 2 2 ...
## $ sem: num [1:384] 1 2 3 4 5 6 7 8 1 2 ...
## $ fil: num [1:384] 1 1 1 1 1 1 1 1 1 1 ...
## $ col: num [1:384] 2 2 2 2 2 2 2 2 3 3 ...
## $ tra: chr [1:384] "A" "A" "A" "A" ...
## $ sil: num [1:384] 0 0 1 1 0 0 1 0 1 3 ...
## $ lab: num [1:384] 0 0 0 0 0 0 6 1 2 0 ...# Transformando filas, columnas y experimentos en factores
da <- transform(da,
fil = factor(fil),
col = factor(col),
exp = factor(exp),
tra = factor(tra))
# Leyendo las seis primeras líneas
head(da)## exp sem fil col tra sil lab
## 1 1 1 1 2 A 0 0
## 2 1 2 1 2 A 0 0
## 3 1 3 1 2 A 1 0
## 4 1 4 1 2 A 1 0
## 5 1 5 1 2 A 0 0
## 6 1 6 1 2 A 0 0
Nota
Qué es una regresión? Y un ANOVA? Cuál es la principal diferencia entre ambos? Qué supuestos estadísticos debemos asumir cuando llevemos a cabo este tipo de análisis? Estas y otras preguntas son críticas en la aplicación de modelos lineales a la resolución de problemas estadísticos.
Esquema conceptual de los pasos que deben seguirse a la hora de ajustar un modelo lineal univariante.(Luis Cayuela, 2015).
Modelo Lineal Generalizado(MLG)
Qué son los MLG?
Los modelos lineales (regresión, ANOVA, ANCOVA), se basan en los siguientes supuestos:
Los errores se distribuyen normalmente.
La varianza es constante.
La variable respuesta se relaciona linealmente con la(s) variable(s) independiente(s).
En muchas ocasiones, sin embargo, nos encontramos con que uno o varios de estos supuestos no se cumplen por la naturaleza de la información.
Estos problemas se pueden llegar a solucionar mediante la transformación de la ariable respuesta (por ejemplo tomando logaritmos).
Sin embargo estas transformaciones no siempre consiguen corregir la falta de normalidad, la heterocedasticidad (varianza no constante) o la no linealidad de nuestros datos.
Además resulta muchas veces interpretar los resultados obtenidos, si utilizamos transformaciones de la variable.
En el Modelo Lineal Generalizado (GLM) tenemos variables de respuesta asociadas a covariables
Mientras en el Modelo Lineal combinamos aditividad de los afectos de las covariables con normalidad de las respuestas y homoscedasticidad, en el GML estas tres cosas no se satisfacen necesariamente.
Los GLM permiten incluir respuestas no normales, como binomial, Poisson, Gamma y otras distibuciones, y en la teoría clásica la estimación se realiza mediante el método de máxima verosimilitud.
Ciertos tipos de variables respuesta sufren invariablemente la violación de estos dos supuestos de los modelos normales y los GLM ofrecen una buena alternativa para tratarlos. Especícamente, podemos considerar utilizar GLM cuando la variable respuesta es:
un conteo de casos (p.e. abundancia de una especie);
un conteo de casos expresados como proporciones (p.e. porcentaje de plántulas muertas en un experimento de vivero);
una respuesta binaria (p.e. vivo o muerto, infectado o no infectado);
Regresando al ejemplo...
#------------------------------------------------------------
# Visualización de los tratamientos durante 8 semanas
# dentro de cada experimento
#------------------------------------------------------------
useOuterStrips(xyplot(lab + sil ~ sem | exp + tra,
data = da,
auto.key = TRUE,
type = c("p", "a")))
db <- aggregate(cbind(lab, sil) ~ exp + tra + col + fil,
data = da, FUN = sum)
head(db)## exp tra col fil lab sil
## 1 3 A 1 1 0 0
## 2 1 B 1 1 0 1
## 3 2 B 1 1 0 0
## 4 1 A 2 1 7 3
## 5 2 C 2 1 10 0
## 6 3 C 2 1 3 0str(db)## 'data.frame': 48 obs. of 6 variables:
## $ exp: Factor w/ 3 levels "1","2","3": 3 1 2 1 2 3 2 1 3 3 ...
## $ tra: Factor w/ 4 levels "A","B","C","D": 1 2 2 1 3 3 1 4 4 2 ...
## $ col: Factor w/ 4 levels "1","2","3","4": 1 1 1 2 2 2 3 3 3 4 ...
## $ fil: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ lab: num 0 0 0 7 10 3 7 0 0 0 ...
## $ sil: num 0 1 0 3 0 0 4 0 0 0 ...Muchos de los métodos estadísticos más comunes, como la t de Student o la regresión, asumen que la varianza es constante, pero en muchas aplicaciones este supuesto no es aplicable. Y es precisamente en estos casos cuando los GLM pueden ser de gran utilidad. Los GLM tienen dos propiedades importantes:
La estructura de los errores.
La función de vínculo.
Regresando al ejemplo...
Antes de ajustar un GLM, vamos a ver cómo sería un modelo de ANOVA clásico ajustado a este problema.
Análisis conjunta de una serie de tres experimentos en diseño cuadrado latino 4X4, considerando la variable insectos silvestres
m_01 <- aov(terms(sil ~ exp/(fil + col) + tra, keep.order = TRUE),
data = db)
par(mfrow = c(2, 2))
plot(m_01)
summary(m_01)## Df Sum Sq Mean Sq F value Pr(>F)
## exp 2 12.87 6.437 1.490 0.246
## exp:fil 9 38.12 4.236 0.980 0.480
## exp:col 9 19.62 2.181 0.505 0.857
## tra 3 27.67 9.222 2.134 0.122
## Residuals 24 103.71 4.321# Test de normalidad de dichos residuos con la función `shapiro.test()`.
shapiro.test(m_01$residuals)##
## Shapiro-Wilk normality test
##
## data: m_01$residuals
## W = 0.9791, p-value = 0.5418El modelo clásico de la ANOVA se ajustó bien a los datos en cuestión
Ahora haremos el análisis conjunta de una serie de tres experimentos en diseño cuadrado latino 4X4, considerando la variable insectos de laboratorio
m_02 <- aov(terms(lab ~ exp/(fil + col) + tra, keep.order = TRUE),
data = db)
par(mfrow = c(2, 2))
plot(m_02)
summary(m_02)## Df Sum Sq Mean Sq F value Pr(>F)
## exp 2 99.5 49.77 2.009 0.156
## exp:fil 9 177.4 19.71 0.795 0.624
## exp:col 9 143.9 15.99 0.645 0.748
## tra 3 66.4 22.14 0.893 0.459
## Residuals 24 594.7 24.78# Test de normalidad de dichos residuos con la función `shapiro.test()`.
shapiro.test(m_02$residuals)##
## Shapiro-Wilk normality test
##
## data: m_02$residuals
## W = 0.94504, p-value = 0.02547El modelo clásico de la ANOVA no se ajustó bien a los datos en cuestión.
Algunas consideraciones ( ejemplo de datos de conteo... ):
Variable respuesta: \(Y_i\) Número de insectos capturados cada semana, durante 8 semanas, en muestras de tamaño \(n=384\).
Distribución: Poisson.
# Visualización de los tratamientos dentro de cada experimento
xyplot(lab + sil ~ tra | exp,
data = da,
auto.key = TRUE,
type = c("p", "a"))
MLG
Supongamos que observamos las variables de respuesta \(Y_1,\ldots,Y_n\) que son v.a. independientes relacionadas con las covariables \(x_{i1},x_{i2},\ldots,x_{ip}\), \(1 \le i \le n\).
Genericamente pensemos en una respuesta \(Y_y\) con covariables \(x_{1},x_{2},\ldots,x_{p}\)
Componentes del Modelo
- Podemos pensar que el GLM posee tres componentes:
- Componente Aleatorio : la variable de respuesta Y tiene distribución
\[ exp\left\{ \frac{y\theta-b(\theta)}{a(\phi)}+c(y,\phi)\right\}, \] siendo \(\theta\) el parámetro canónico, \(\phi\) es un parámetro de dispersión y las funciones \(a()\), \(b()\), y \(c()\) son conocidas. Además, se cumple que:
\(\mu=E(Y)=b'(\theta)\) y \(Var(Y)=a(\phi)b''(\theta)\)
Componente Sistemático : el vector de covariables \(\mathbf{x'}=(x_1,x_2,\ldots,x_p)\) que da origen al predictor lineal \[ \eta=\sum_{j=1}^px_j\beta_j=\mathbf{x'\beta} \] siendo \(\mathbf{\beta}\) el vector a estimar
Función de enlace o link : relaciona los dos componentes \(\mu\) y \(\eta\) \[ g(\mu)=\eta \]
Nota: En algunos casos \(a(\phi)\) es de la forma \(a(\phi)=\frac{\phi}{w}\), siendo w un peso conocido.
Los modelos lineales generalizados permiten dos extensiones:
I. podemos tratar distribuiciones que pertenezcan a una familia exponencial
- podemos elegir una función de enlace que sea una función monótona y diferenciable.
El Modelo Lineal Generalizado tuvo mucha difusión a partir del libro de McCullagh y Nelder (1989). En estos modelos la variable de respuesta \(Y_i\) sigue una distribución que pertenece a una familia exponencial con media \(\mu_i\) que es una función, por lo general no lineal, de \(\mathbf{x_i'\beta}\).
Nota: Es posible obtener el modelo clásico de ANOVA utilizando la función siguiente:
#############################################
# ANOVA para a variable insectos silvestres #
#############################################
m0_2 <- glm(terms(sil ~ exp/(fil + col) + tra, keep.order = TRUE),
data = db, family = gaussian)
#summary(m0_2)
1-pchisq(103.71/4.321181,24)## [1] 0.4615754anova(m0_2, test = "F")## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: sil
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 47 202.00
## exp 2 12.875 45 189.12 1.4898 0.2455
## exp:fil 9 38.125 36 151.00 0.9803 0.4804
## exp:col 9 19.625 27 131.38 0.5046 0.8566
## tra 3 27.667 24 103.71 2.1342 0.1223shapiro.test(m0_2$residuals)##
## Shapiro-Wilk normality test
##
## data: m0_2$residuals
## W = 0.9791, p-value = 0.5418####################################################
## https://www.ime.usp.br/~giapaula/
##
####################################################
#------------------------------------------------------------#
par(mfrow=c(1,1))
X <- model.matrix(m0_2)
n <- nrow(X)
p <- ncol(X)
H <- X%*%solve(t(X)%*%X)%*%t(X)
h <- diag(H)
si <- lm.influence(m0_2)$sigma
r <- resid(m0_2)
tsi <- r/(si*sqrt(1-h))
#
ident <- diag(n)
epsilon <- matrix(0,n,100)
e <- matrix(0,n,100)
e1 <- numeric(n)
e2 <- numeric(n)
#
for(i in 1:100){
epsilon[,i] <- rnorm(n,0,1)
e[,i] <- (ident - H)%*%epsilon[,i]
u <- diag(ident - H)
e[,i] <- e[,i]/sqrt(u)
e[,i] <- sort(e[,i]) }
#
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(tsi,e1,e2)
#
par(pty="s")
qqnorm(tsi,xlab="Percentil da N(0,1)",
ylab="Residuo Studentizado", ylim=faixa, pch=16, main="")
par(new=TRUE)
qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1, main="")
par(new=TRUE)
qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1, main="")
par(new=TRUE)
qqnorm(med,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=2, main="")
#------------------------------------------------------------#hnp(m0_2, print=TRUE,verb=TRUE,halfnormal = T,xlab="Cuantiles teóricos")
#################################################
# ANOVA para a variable insectos de laboratorio #
#################################################
m0_3 <- glm(terms(lab ~ exp/(fil + col) + tra, keep.order = TRUE),
data = db, family = gaussian)
#summary(m0_3)
1-pchisq(594.71/24.77951,24)## [1] 0.4615933anova(m0_3, test = "F")## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: lab
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 47 1081.92
## exp 2 99.542 45 982.38 2.0085 0.1561
## exp:fil 9 177.375 36 805.00 0.7953 0.6238
## exp:col 9 143.875 27 661.12 0.6451 0.7480
## tra 3 66.417 24 594.71 0.8934 0.4588####################################################
## https://www.ime.usp.br/~giapaula/
##
####################################################
#------------------------------------------------------------#
par(mfrow=c(1,1))
X <- model.matrix(m0_3)
n <- nrow(X)
p <- ncol(X)
H <- X%*%solve(t(X)%*%X)%*%t(X)
h <- diag(H)
si <- lm.influence(m0_3)$sigma
r <- resid(m0_3)
tsi <- r/(si*sqrt(1-h))
#
ident <- diag(n)
epsilon <- matrix(0,n,100)
e <- matrix(0,n,100)
e1 <- numeric(n)
e2 <- numeric(n)
#
for(i in 1:100){
epsilon[,i] <- rnorm(n,0,1)
e[,i] <- (ident - H)%*%epsilon[,i]
u <- diag(ident - H)
e[,i] <- e[,i]/sqrt(u)
e[,i] <- sort(e[,i]) }
#
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(tsi,e1,e2)
#
par(pty="s")
qqnorm(tsi,xlab="Percentil da N(0,1)",
ylab="Residuo Studentizado", ylim=faixa, pch=16, main="")
par(new=TRUE)
qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1, main="")
par(new=TRUE)
qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1, main="")
par(new=TRUE)
qqnorm(med,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=2, main="")
#------------------------------------------------------------#hnp(m0_3, print=TRUE,verb=TRUE,halfnormal = T,xlab="Cuantiles teóricos")
Regresando al ejemplo...
#---------------------------------------------------------------#
# Análisis de la suma para el experimento 1, considerando la variable #
#insectos silvestres #
#---------------------------------------------------------------#
m1 <- glm(sil ~ fil + col + tra,
data = subset(db, exp == "1"),
family = gaussian)
#summary(m1)
1-pchisq(14.50/2.416667,6)## [1] 0.4231902shapiro.test(m1$residuals)##
## Shapiro-Wilk normality test
##
## data: m1$residuals
## W = 0.96489, p-value = 0.7507anova(m1, test = "F")## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: sil
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 15 45.75
## fil 3 4.25 12 41.50 0.5862 0.6459
## col 3 6.25 9 35.25 0.8621 0.5101
## tra 3 20.75 6 14.50 2.8621 0.1264
hnp(m1, print=TRUE,verb=TRUE,halfnormal = T,xlab="Cuantiles teóricos")
#---------------------------------------------------------------------#
# Análisis de la suma para el experimento 2, considerando la variable #
#insectos silvestres #
#---------------------------------------------------------------------#
m2 <- glm(sil ~ fil + col + tra,
data = subset(db, exp == "2"),
family = gaussian)
#summary(m2)
1-pchisq(40.875/6.8125,6)## [1] 0.4231901shapiro.test(m2$residuals)##
## Shapiro-Wilk normality test
##
## data: m2$residuals
## W = 0.95161, p-value = 0.5155anova(m2, test = "F")## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: sil
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 15 104.438
## fil 3 19.688 12 84.750 0.9633 0.4687
## col 3 8.688 9 76.062 0.4251 0.7422
## tra 3 35.188 6 40.875 1.7217 0.2613
hnp(m2, print=TRUE,verb=TRUE,halfnormal = F,xlab="Cuantiles teóricos")
#---------------------------------------------------------------------#
# Análisis de la suma para el experimento 2, considerando la variable #
#insectos silvestres #
#---------------------------------------------------------------------#
m3 <- glm(sil ~ fil + col + tra,
data = subset(db, exp == "3"),
family = gaussian)
#summary(m3)
1-pchisq(17.875/2.979167,6)## [1] 0.4231902shapiro.test(m3$residuals)##
## Shapiro-Wilk normality test
##
## data: m3$residuals
## W = 0.9627, p-value = 0.711anova(m3, test = "F")## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: sil
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 15 38.938
## fil 3 14.1875 12 24.750 1.5874 0.2880
## col 3 4.6875 9 20.062 0.5245 0.6813
## tra 3 2.1875 6 17.875 0.2448 0.8623
hnp(m3, print=TRUE,verb=TRUE,halfnormal = F,xlab="Cuantiles teóricos")
#---------------------------------------------------------------------#
# Análisis de la suma para el experimento 1, considerando la variable #
#insectos de laboratorio #
#---------------------------------------------------------------------#
m1.1 <- glm(lab ~ fil + col + tra,
data = subset(db, exp == "1"),
family = gaussian)
#summary(m1.1)
1-pchisq(33/5.5,6)## [1] 0.4231901shapiro.test(m1.1$residuals)##
## Shapiro-Wilk normality test
##
## data: m1.1$residuals
## W = 0.91736, p-value = 0.153anova(m1.1, test = "F")## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: lab
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 15 114.0
## fil 3 3.5 12 110.5 0.2121 0.8845
## col 3 37.5 9 73.0 2.2727 0.1803
## tra 3 40.0 6 33.0 2.4242 0.1639
hnp(m1.1, print=TRUE,verb=TRUE,halfnormal = F,xlab="Cuantiles teóricos")
#---------------------------------------------------------------------#
# Análisis de la suma para el experimento 2, considerando la variable #
#insectos de laboratorio #
#---------------------------------------------------------------------#
m2.1 <- glm(lab ~ fil + col + tra,
data = subset(db, exp == "2"),
family = gaussian)
#summary(m2.1)
1-pchisq(414.88/69.14583,6)## [1] 0.4231819shapiro.test(m2.1$residuals)##
## Shapiro-Wilk normality test
##
## data: m2.1$residuals
## W = 0.91509, p-value = 0.1406anova(m2.1, test = "F")## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: lab
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 15 781.44
## fil 3 143.187 12 638.25 0.6903 0.5904
## col 3 88.188 9 550.06 0.4251 0.7422
## tra 3 135.188 6 414.88 0.6517 0.6104
hnp(m2.1, print=TRUE,verb=TRUE,halfnormal = F,xlab="Cuantiles teóricos")
#---------------------------------------------------------------------#
# Análisis de la suma para el experimento 3, considerando la variable #
#insectos de laboratorio #
#---------------------------------------------------------------------#
m3.1 <- glm(lab ~ fil + col + tra,
data = subset(db, exp == "3"),
family = gaussian)
#summary(m3.1)
1-pchisq(36.875/6.145833,6)## [1] 0.42319anova(m3.1, test = "F")## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: lab
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 15 86.938
## fil 3 30.6875 12 56.250 1.6644 0.2723
## col 3 18.1875 9 38.062 0.9864 0.4598
## tra 3 1.1875 6 36.875 0.0644 0.9768
hnp(m3.1, print=TRUE,verb=TRUE,halfnormal = F,xlab="Cuantiles teóricos")
Modelos de cuasi-verosimilitud (quasi)
En todas las familias, la varianza de la respuesta dependerá de la media, y el parámetro de escala actuará como un multiplicador. La forma de dependencia de la varianza respecto de la media es una característica de la distribución de respuesta, por ejemplo para la distribución de Poisson será \(Var[Y] = \mu\).
Para estimación e inferencia de cuasi-verosimilitud, la distribución de respuesta precisa no está especificada, sino más bien sólo la función de enlace y la forma en que la función de varianza depende de la media. Puesto que la estimación de cuasi-verosimilitud utiliza formalmente las mismas técnicas de la distribución gaussiana, esta familia permite ajustar modelos gaussianos con funciones de enlace no estándar o incluso con funciones de varianza.
Función de Varianza para el modelo de Poisson
Este modelo asume que \[ E[Y_i]=Var[Y_i]=\mu_i \] sin embargo es posible que un conjunto de datos tengan una dispersión mayor. Cuando los datos exhiben sobredispersión, se puede tomar uno de los siguientes caminos:
Suponer que \(Var[Y_i]=\sigma^2\mu_i\) y estimar \(\sigma^2\) usando un modelo de quasi- verosimilitud, como en el caso binonial.
Sumergir a la variable de respuesta en una familia de distribuciones que contemple una dispersión mayor: Binomial Negativa
Ahora haremos la análisis conjunta de una serie de tres experimentos en diseño cuadrado latino 4X4, considerando la variable insectos silvestres y de laboratorio
#######################
# Insectos silvestres #
#######################
attach(db)
m0.1 <- aov(terms(sil ~ exp/(fil + col) + tra, keep.order = TRUE),
data = db)
summary(m0.1)## Df Sum Sq Mean Sq F value Pr(>F)
## exp 2 12.87 6.437 1.490 0.246
## exp:fil 9 38.12 4.236 0.980 0.480
## exp:col 9 19.62 2.181 0.505 0.857
## tra 3 27.67 9.222 2.134 0.122
## Residuals 24 103.71 4.321#Calidad del ajuste basado en shapiro wilk
shapiro.test(m0.1$residuals)##
## Shapiro-Wilk normality test
##
## data: m0.1$residuals
## W = 0.9791, p-value = 0.5418attach(db)## The following objects are masked from db (pos = 3):
##
## col, exp, fil, lab, sil, tra# Teste de Tukey p
dados.anova1<-anova(m0.1)
tuk1<-HSD.test(sil,tra,dados.anova1$Df[5],
dados.anova1$Mean[5],alpha=0.05)
tuk1## $statistics
## MSerror Df Mean CV MSD
## 4.321181 24 1.5 138.583 2.341077
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey tra 4 3.901262 0.05
##
## $means
## sil std r Min Max Q25 Q50 Q75
## A 2.750000 2.701010 12 0 9 0.75 3 3.25
## B 1.083333 1.831955 12 0 6 0.00 0 1.25
## C 0.750000 1.215431 12 0 4 0.00 0 1.00
## D 1.416667 1.928652 12 0 6 0.00 1 1.50
##
## $comparison
## NULL
##
## $groups
## sil groups
## A 2.750000 a
## D 1.416667 a
## B 1.083333 a
## C 0.750000 a
##
## attr(,"class")
## [1] "group"
La función glht funciona con Hipótesis lineales generales y comparaciones múltiples para modelos paramétricos, incluidos modelos lineales generalizados, modelos de efectos mixtos lineales y modelos de supervivencia.
mult.comp.m0.1 <- glht(m0.1, linfct=mcp(tra="Tukey"))
summary(mult.comp.m0.1)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = terms(sil ~ exp/(fil + col) + tra, keep.order = TRUE),
## data = db)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 -1.6667 0.8486 -1.964 0.229
## C - A == 0 -2.0000 0.8486 -2.357 0.113
## D - A == 0 -1.3333 0.8486 -1.571 0.413
## C - B == 0 -0.3333 0.8486 -0.393 0.979
## D - B == 0 0.3333 0.8486 0.393 0.979
## D - C == 0 0.6667 0.8486 0.786 0.860
## (Adjusted p values reported -- single-step method)
###########################
# Insectos de laboratorio #
###########################
m0.2 <- aov(terms(lab ~ exp/(fil + col) + tra, keep.order = TRUE),
data = db)
#summary(m0.2)
#El modelo está bien ajustado
#Calidad del ajuste basado en shapiro wilk
shapiro.test(m0.1$residuals)##
## Shapiro-Wilk normality test
##
## data: m0.1$residuals
## W = 0.9791, p-value = 0.5418attach(db)## The following objects are masked from db (pos = 3):
##
## col, exp, fil, lab, sil, tra## The following objects are masked from db (pos = 4):
##
## col, exp, fil, lab, sil, tra# Teste de Tukey p
dados.anova2<-anova(m0.2)
tuk2<-HSD.test(lab,tra,dados.anova2$Df[5],
dados.anova2$Mean[5],alpha=0.05)
tuk2## $statistics
## MSerror Df Mean CV MSD
## 24.77951 24 3.291667 151.2274 5.606101
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey tra 4 3.901262 0.05
##
## $means
## lab std r Min Max Q25 Q50 Q75
## A 4.666667 7.900901 12 0 28 0 1.5 6.25
## B 1.416667 1.880925 12 0 5 0 0.5 2.50
## C 3.583333 3.203928 12 0 10 1 2.5 5.00
## D 3.500000 4.011348 12 0 12 0 2.0 5.50
##
## $comparison
## NULL
##
## $groups
## lab groups
## A 4.666667 a
## C 3.583333 a
## D 3.500000 a
## B 1.416667 a
##
## attr(,"class")
## [1] "group"
mult.comp.m0.2 <- glht(m0.2, linfct=mcp(tra="Tukey"))
summary(mult.comp.m0.2)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = terms(lab ~ exp/(fil + col) + tra, keep.order = TRUE),
## data = db)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 -3.25000 2.03222 -1.599 0.398
## C - A == 0 -1.08333 2.03222 -0.533 0.950
## D - A == 0 -1.16667 2.03222 -0.574 0.939
## C - B == 0 2.16667 2.03222 1.066 0.713
## D - B == 0 2.08333 2.03222 1.025 0.737
## D - C == 0 -0.08333 2.03222 -0.041 1.000
## (Adjusted p values reported -- single-step method)
Ejemplo 2
Ejemplo 7.1, página 121, del libro: Modelos de Regresión Gauss M. Cordero y Clarice G. B. Demetrio
Las muestras de 20 insectos Heliothis virescens (plaga de algodón) resistentes a la cipermetrina, se expusieron a dosis crecientes del insecticida, dos días después de la emergencia de la pupa (Collet, 2002). Después de 72 horas se contó el número de insectos muertos.
Cyper = Cyper = read.table("Insectos.txt",header=TRUE)
head(Cyper)## sex dose ldose y
## 1 M 1 0 1
## 2 M 2 1 4
## 3 M 4 2 9
## 4 M 8 3 13
## 5 M 16 4 18
## 6 M 32 5 20Reconoce como un objeto aislado cada columna de un Data.frame
attach(Cyper)Antes de ajustar un GLM, vamos a ver cómo sería un modelo lineal normal ajustado a este problema.
fit1 = lm(y~dose+sex)
summary(fit1)##
## Call:
## lm(formula = y ~ dose + sex)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9483 -2.5275 0.1047 2.1198 4.3385
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.26741 1.66750 1.360 0.206989
## dose 0.51421 0.08981 5.726 0.000285 ***
## sexM 3.16667 1.94491 1.628 0.137928
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.369 on 9 degrees of freedom
## Multiple R-squared: 0.7975, Adjusted R-squared: 0.7524
## F-statistic: 17.72 on 2 and 9 DF, p-value: 0.0007574fit2 = lm(y~dose)
summary(fit2)##
## Call:
## lm(formula = y ~ dose)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.3650 -3.0006 -0.1919 2.2998 5.9218
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.85075 1.46213 2.634 0.025005 *
## dose 0.51421 0.09694 5.305 0.000345 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.636 on 10 degrees of freedom
## Multiple R-squared: 0.7378, Adjusted R-squared: 0.7116
## F-statistic: 28.14 on 1 and 10 DF, p-value: 0.0003453shapiro.test(fit2$residuals)##
## Shapiro-Wilk normality test
##
## data: fit2$residuals
## W = 0.93692, p-value = 0.4592Algunas consideraciones ( ejemplo del Insectos... ):
Variable respuesta: \(Y_i\) Número de insectos muertos en muestras de tamaño \(m_i=20\)
Distribución: Binomial
En estadística, la distribución binomial o distribución binómica es una distribución de probabilidad discreta que cuenta el número de éxitos en una secuencia de n ensayos de Bernoulli independientes entre sí, con una probabilidad fija p de ocurrencia del éxito entre los ensayos.
plot(ldose, y/20, pch=c(rep("*",6), rep("+",6)), col=c(rep("green",6),
rep("red",6)), xlab="log(dosis)", ylab="Proporción de insectos muertos")
legend("bottomright", c("Macho", "Hembra"), pch=c("*","+"),
col=c("green","red"))
title("Proporciones observadas y las curvas ajustadas")
Se observa que las hembras son más resistentes, mientras que para matar al 100% de las hembras hay necesidad de una dosis 2 veces mayor que para matar al 100% de los machos. Puede verificarse que la dosis letal para p = 0.9 para las hembras está fuera del rango estudiado, lo cual es peligroso, mientras que para la dosis mayor que 32 no se sabe si el comportamiento es el mismo.
Regresando al ejemplo...
Si el interés se centra en la estimación de esa dosis ès necesario aumentar la amplitud de la dosis para las hembras en un nuevo experimento.
resp<-cbind(y,20-y);resp## y
## [1,] 1 19
## [2,] 4 16
## [3,] 9 11
## [4,] 13 7
## [5,] 18 2
## [6,] 20 0
## [7,] 0 20
## [8,] 2 18
## [9,] 6 14
## [10,] 10 10
## [11,] 12 8
## [12,] 16 4La comparación de algunos modelos GLM considerando la distribución binomial
mod1 = glm(resp~1, family=binomial)
summary(mod1)##
## Call:
## glm(formula = resp ~ 1, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.9833 -2.7215 0.1118 2.0401 5.5538
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1503 0.1295 -1.161 0.246
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 124.88 on 11 degrees of freedom
## Residual deviance: 124.88 on 11 degrees of freedom
## AIC: 156.99
##
## Number of Fisher Scoring iterations: 3mod2 = glm(resp~ldose, family=binomial)
summary(mod2)##
## Call:
## glm(formula = resp ~ ldose, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7989 -0.8267 -0.1871 0.8950 1.9850
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7661 0.3701 -7.473 7.82e-14 ***
## ldose 1.0068 0.1236 8.147 3.74e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 124.876 on 11 degrees of freedom
## Residual deviance: 16.984 on 10 degrees of freedom
## AIC: 51.094
##
## Number of Fisher Scoring iterations: 4mod3 = glm(resp~sex, family=binomial)
summary(mod3)##
## Call:
## glm(formula = resp ~ sex, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.7887 -2.9371 0.1015 2.3400 4.9522
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4754 0.1878 -2.532 0.0113 *
## sexM 0.6425 0.2623 2.449 0.0143 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 124.88 on 11 degrees of freedom
## Residual deviance: 118.80 on 10 degrees of freedom
## AIC: 152.91
##
## Number of Fisher Scoring iterations: 4mod4 = glm(resp~ldose+sex, family=binomial)
summary(mod4)##
## Call:
## glm(formula = resp ~ ldose + sex, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.10540 -0.65343 -0.02225 0.48471 1.42944
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.4732 0.4685 -7.413 1.23e-13 ***
## ldose 1.0642 0.1311 8.119 4.70e-16 ***
## sexM 1.1007 0.3558 3.093 0.00198 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 124.8756 on 11 degrees of freedom
## Residual deviance: 6.7571 on 9 degrees of freedom
## AIC: 42.867
##
## Number of Fisher Scoring iterations: 4Comparación de modelos por el Criterio de información de Akaike (AIC)
AIC(mod1,mod2,mod3,mod4)## df AIC
## mod1 1 156.98599
## mod2 2 51.09443
## mod3 2 152.90901
## mod4 3 42.86747
Cuando se tiene una serie de modelos \(M_1\), \(M_2\), ... con parámetros \(\theta_1\), \(\theta_2\), ... respectivamente, una metodología para compararlos corresponde a la función de máxima verosimilitud (likelihood). La máxima verosimilitud permite seleccionar el modelo que realiza el mejor ajuste de los datos pero no penaliza su complejidad, lo que si sucede cuando se emplean medidas de contraste comoel AIC.
EL criterio hace uso del Log-likelihood (log Lik), que es el logaritmo de máxima verosimilitud, y sustraen un término proporcional al número de parámetros \((\theta)\) en el modelo
\[ AIC=-2\times log Lik + 2\theta \]
Los mejores modelos son aquellos que presentaron el menor valor de AIC
Ejemplos...
- Normal: \(Y \sim N(\mu,\sigma)\) \[ f(y,\theta,\phi)=\frac{1}{\sqrt{2\pi\sigma^2}}exp\left(-\frac{1}{2}\frac{(y-\mu)^2}{\sigma^2} \right)\\ =exp\left( \frac{y\mu-\mu^2/2}{\sigma^2}-\frac{1}{2}\left[ \frac{y^2}{\sigma^2}+log(2\pi\sigma^2)\right] \right), \] por lo tanto \(\theta=\mu\), \(b(\theta)=\frac{\mu^2}{2}\), \(\phi=\sigma^2\), \(a(\phi)=\phi\) y \(c(y,\phi)=-\frac{1}{2}\left[\frac{y^2}{\phi}+log(2\pi\phi) \right]\).
\[ E(Y)=\mu \]
En el caso heteroscedástico \(Y\sim N(\mu,\frac{\sigma^2}{w})\), donde w es un peso conocido, tenemos \(\phi=\sigma^2\) y \(a(\phi)=\frac{\phi}{w}\).
- Caso Binomial: \(Y\sim Bi(n,p)\)
Considerando \(\frac{Y}{n}=\) , la proporción de éxitos.
\[ P\left(\frac{Y}{n}=y \right)=P(Y=ny)=\binom{n}{ny}p^{ny}(1-p)^{n-ny} \] \[ P\left(\frac{Y}{n}=y \right)=exp\left(\frac{ylog\left(\frac{p}{1-p} \right)+log(1-p)}{1/n} +log\binom{n}{ny}\right), \] por lo tanto \(\theta=log\frac{p}{1-p}\) , \(b(\theta)=log(1+e^{\theta})\) , \(\phi=n\) , \(a(\phi)=\frac{1}{\phi}\) y \(c(y,\phi)=\binom{\phi}{\phi y}\) .
\[ E\left(\frac{Y}{n} \right)=p=\frac{e^\theta}{1+e^\theta}=\frac{1}{1+e^{-\theta}} \]
- Caso Poisson: \(Y\sim P(\lambda)\) \[ P(Y=y)=e^{-\lambda}\frac{\lambda^y}{y!}\\ P(Y=y)=exp(ylog\lambda-\lambda-log(y!)) \] por lo tanto \(\theta=log\lambda\), \(b(\theta)=e^{\theta}\), \(\phi=1\) y \(c(y,\phi)=-log(y!)\) \[ E(Y)=\lambda=e^{\theta} \]
Función de enlace o link:
Esta función relaciona el predictor lineal \(\eta\) con la esperanza \(\mu\) de la respuesta Y. A diferencia del modelo lineal clásico, aquí introducimos una función uno-a-uno continua y diferenciable, \(g(\mu)\), tal que \[ \eta=g(\mu) \] En el caso Binomial, por ejemplo, tenemos que \(\mu \in (0,1)\) y la función de enlace tiene que pertenecer a la recta de los números reales. Suelen usarse 3 funciones de enlaces:
- Logit: \(\eta=log\frac{\mu}{1-\mu}\)
- Probit: \(\eta=\Phi^{-1}(\mu)\)
- Complemento log-log: \(\eta=log(-log(1-\mu))\)
Links Canónicos:
En el caso normal mostramos que si \(Y\sim(\mu,\sigma^2)\) el parámetro canónico es \(\theta=\mu\).
En el caso binomial \(Y\sim Bi(n,p)\) en el que consideramos \(\frac{Y}{n}\) vimos que el vector canónico es \(\theta=logit(\pi)\). Estos son los enlaces más usados en cada caso.
Cuando usamos \(\eta=\theta\) el modelo tiene el enlace canónico convencional. Es conveniente usar el enlace convencional, ya que algunas cosas se simplifican, pero la posibilidad de usarlo dependerá de los datos con los que estemos trabajando.
- Normal: \(\eta\mu\)
- Poisson: \(\eta=log\mu\)
- Binomial: \(\eta=log\frac{\mu}{1-\mu}\)
- Gamma: \(\eta=\mu^{-1}\)
Comprobamos la calidad del ajuste basandonos en la deviance.
Modelo1
m1 = 1-pchisq(mod1$deviance,mod1$df.residual); m1## [1] 0Modelo2
m2 = 1-pchisq(mod2$deviance,mod2$df.residual); m2## [1] 0.07471804Modelo3
m3 = 1-pchisq(mod3$deviance,mod3$df.residual); m3## [1] 0Modelo4
m4 = 1-pchisq(mod4$deviance,mod4$df.residual); m4## [1] 0.6623957Hay evidencia significativa de que el ajuste del modelo 4 és apropiado (p=0.66)
Función de Verosimilitud para el GLM:
Sea Y una v.a. con función de densidad de probabilidad perteneciente a una familia exponencial dada por: \[ f_Y(y,\theta,\phi)=exp\left\{ \frac{y\theta-b(\theta)}{a(\phi)}+c(y,\phi)\right\}, \] para algunas funciones conocidas \(a(\phi)\), \(b(\theta)\) y \(c(y , \phi)\). \(S_i\) es un parámetro conocido, esta es una familia exponencial con parámetro canónico natural \(\theta\).
Si \(\phi\) no es conocido, ésta puede ser una familia exponencial en \((\theta, \phi)\) o no. \(\phi\) es un parámetro de dispersión o de forma.
La media \(E(Y)\) es sólo función de \(\theta\) y es por lo tanto el parámetro de interés; \(\phi\) en general es tratado como un parámetro nuisance o de ruido. En la mayoría de los casos \(\phi\) no será tratado tal como es tratado \(\theta\). Estimaremos y haremos inferencia bajo un valor asumido de \(\phi\) y si \(\phi\) necesita ser estimado, lo estimaremos y luego será tomado como un valor fijo y conocido.
Esta familia incluye distribuciones simétricas, asimétricas, discretas y continuas, tales como la distribución Normal, Binomial, Poisson o Gamma.
Momentos de una familia exponencial:
Deduciremos el primer y segundo momento de una familia exponencial a partir del logaritmo de su verosimilitud. \[ L(\theta,y)=\frac{y\theta-b(\theta)}{a(\phi)}+c(y,\phi). \] Su primera derivada es: \[ L'(\theta,y)=\frac{\partial L(\theta,y)}{\partial\theta}=\frac{y-b'(\theta)}{a(\phi)}, \] mientras que su derivada segunda es: \[ L''(\theta,y)=\frac{\partial L^2(\theta,y)}{\partial^2\theta}=\frac{-b''(\theta)}{a(\phi)}. \] Como \(E\left( \frac{\partial L(\theta,y)}{\partial\theta}\right)=0\), entonces \[ 0=E(L'(\theta,y))=E\left[ \frac{y-b'(\theta)}{a(\phi)}\right] \] y por lo tanto \[ \mu=E(Y)=b'(\theta). \] Además, sabemos que \[ E(L''(\theta,y))=-E[(L'(\theta,y))^2], \] entonces \[ Var(L'(\theta,y))=E[(L'(\theta,y))^2]=-E(L''(\theta,y))=\frac{b''(\theta)}{a(\phi)}. \] Por otro lado, \[ Var(L'(\theta,y))=Var\left(\frac{y-b'(\theta)}{a(\phi)}\right)=\frac{1}{a^2(\phi)}Var(Y) \] y en consecuencia \[ Var(Y)=a(\phi)b''(\theta). \] La varianza es el producto de dos funciones: una que depende del parámetro natural, \(\theta\), y otra que depende sálo del parámetro nuisance \(\phi\). \(V(\theta) = b''(\theta)\) es llamada función de varianza del modelo.
Resumiendo: \[ E(Y)=b'(\theta)\\ Var(Y)=a(\phi)b''(\theta) \] Nota: Estimación de los parámetros: Método de Newton-Raphson y Fisher-scoring
Regresando al ejemplo...:
Vamos a seguir el diagnóstico y comprobar si hay valores atípicos y puntos influyentes.
influenceIndexPlot(mod4,vars=c('Studentized','Cook','Hat'),id.n=3)
plot(c(1,32), c(0,1), type="n", xlab="log(dose)",
ylab="Proporción de insectos muertos", log="x")
points(2**ldose, y/20, pch=c(rep("*",6),rep("+",6)),
col=c(rep("green",6),rep("red",6)))
ld<-seq(0,5,0.1)
lines(2**ld, predict(mod4,data.frame(ldose=ld,
sex=factor(rep("M",length(ld)),levels=levels(sex))),
type="response"), col="green")
lines(2**ld, predict(mod4,data.frame(ldose=ld,
sex=factor(rep("F",length(ld)),levels=levels(sex))),
type="response"), col="red")
title("Proporciones observadas y modelo ajustado")
legend("bottomright", c("Macho", "Hembra"), pch=c("*","+"),
col=c("green","red"))
#https://www.ime.usp.br/~giapaula/envel_bino
#------------------------------------------------------------#
par(mfrow=c(1,1))
X <- model.matrix(mod4)
n <- nrow(X)
p <- ncol(X)
w <- mod4$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
td <- resid(mod4,type="deviance")/sqrt(1-h)
e <- matrix(0,n,100)
#
for(i in 1:100){
dif <- runif(n) - fitted(mod4)
dif[dif >= 0 ] <- 0
dif[dif<0] <- 1
nresp <- dif
fit <- glm(nresp ~ X, family=binomial)
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 de la N(0,1)",
ylab="Componente de la Desviacion", ylim=faixa, pch=16, main="")
#
par(new=TRUE)
#
qqnorm(e1,axes=FALSE,xlab="",ylab="",type="l",ylim=faixa,lty=1, main="")
par(new=TRUE)
qqnorm(e2,axes=FALSE,xlab="",ylab="", type="l",ylim=faixa,lty=1, main="")
par(new=TRUE)
qqnorm(med,axes=FALSE,xlab="", ylab="", type="l",ylim=faixa,lty=2, main="")
Ejemplo: ##Número de insectos por parcela
Universidad de San Carlos de Guatemala
Centro de Telemática (CETE)
Facultad de Agronomía USAC
Use R!
Ejemplo 3
El escarabajo del pino de montaña (Dendroctonus ponderosae) es una especie que posee un exoesqueleto duro y negro y mide unos 5 milímetros. Es una plaga importante de los bosques de pino ponderosa en el estado de Colorado. Los brotes infecciosos de estos escarabajos matan cantidades considerables de pinos.
La extracción de todos los objetos
rm(list=ls(all=TRUE))Comienza cargando las librerías y preparando los datos a utilizar.
datos <- read_excel("Datos_conteo.xlsx",
sheet = "Data_R4")
head(datos)## # A tibble: 6 x 6
## y x1 x2 x3 x4 x5
## <dbl> <chr> <chr> <chr> <chr> <chr>
## 1 0 40.5 37.0 7.50 1.38 1.4
## 2 87 37.0 34.0 6.50 1.72 34.0
## 3 42 57.5 52.5 12.60 0.83 6.2
## 4 4 32.5 32.0 4.00 1.61 16.6
## 5 0 37.5 34.5 7.00 1.78 1.0
## 6 40 44.0 39.8 11.00 1.68 12.2dim(datos)## [1] 44 6str(datos)## tibble [44 x 6] (S3: tbl_df/tbl/data.frame)
## $ y : num [1:44] 0 87 42 4 0 40 0 0 0 45 ...
## $ x1: chr [1:44] "40.5" "37.0" "57.5" "32.5" ...
## $ x2: chr [1:44] "37.0" "34.0" "52.5" "32.0" ...
## $ x3: chr [1:44] "7.50" "6.50" "12.60" "4.00" ...
## $ x4: chr [1:44] "1.38" "1.72" "0.83" "1.61" ...
## $ x5: chr [1:44] "1.4" "34.0" "6.2" "16.6" ...attach(datos) ## The following object is masked from Cyper:
##
## ydatos$y=as.numeric(datos$y)
datos$x1=as.numeric(datos$x1)
datos$x2=as.numeric(datos$x2)
datos$x3=as.numeric(datos$x3)
datos$x4=as.numeric(datos$x4)
datos$x5=as.numeric(datos$x5)- Variable respuesta: \(Y_i\) Número de insectos por parcela (datos de conteo)
-Distribución: Poisson
# Ajuste considerando la distribución de Poisson
fit1 <- glm(y~x1+x2+x3+x4+x5, family=poisson,data=datos)
summary(fit1) ##
## Call:
## glm(formula = y ~ x1 + x2 + x3 + x4 + x5, family = poisson, data = datos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -9.360 -3.991 -2.825 2.102 8.715
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.561426 1.452700 -7.270 3.59e-13 ***
## x1 0.808590 0.050237 16.096 < 2e-16 ***
## x2 -0.470453 0.025382 -18.535 < 2e-16 ***
## x3 -0.726792 0.073200 -9.929 < 2e-16 ***
## x4 1.856551 0.351263 5.285 1.25e-07 ***
## x5 0.071627 0.002985 23.997 < 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: 2377.4 on 43 degrees of freedom
## Residual deviance: 1134.2 on 38 degrees of freedom
## AIC: 1281.1
##
## Number of Fisher Scoring iterations: 7#Comprobamos la calidad del ajuste basado en la deviance.
ajuste1 = 1-pchisq(fit1$deviance,fit1$df.residual); ajuste1## [1] 0#Hay evidencia significativa que el ajuste del modelo no es apropiado (p<0.0001)Vamos a ver el sobre de simulación
#https://www.ime.usp.br/~giapaula/envel_pois
par(mfrow=c(1,1))
X <- model.matrix(fit1)
n <- nrow(X)
p <- ncol(X)
w <- fit1$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
td <- resid(fit1,type="deviance")/sqrt((1-h))
e <- matrix(0,n,100)
#
for(i in 1:100){
nresp <- rpois(n, fitted(fit1))
fit <- glm(nresp ~ X, family=poisson)
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 de la N(0,1)",
ylab="Componente de la Desviacion", ylim=faixa, pch=16, main="")
par(new=TRUE)
#
qqnorm(e1,axes=FALSE,xlab="",ylab="",type="l",ylim=faixa,lty=1, main="")
par(new=TRUE)
qqnorm(e2,axes=FALSE,xlab="",ylab="", type="l",ylim=faixa,lty=1, main="")
par(new=TRUE)
qqnorm(med,axes=FALSE,xlab="", ylab="", type="l",ylim=faixa,lty=2, main="")
#------------------------------------------------------------# Veamos el problema de sobredispersión
fit2 <- glm.nb(y~x1+x2+x3+x4+x5, data=datos)
summary(fit2)##
## Call:
## glm.nb(formula = y ~ x1 + x2 + x3 + x4 + x5, data = datos, init.theta = 0.3263879211,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.99661 -1.10980 -0.67357 -0.05588 1.47701
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.93258 12.73411 -1.644 0.1002
## x1 1.00489 0.41770 2.406 0.0161 *
## x2 -0.38131 0.24540 -1.554 0.1202
## x3 -0.94488 0.61096 -1.547 0.1220
## x4 1.49305 3.07811 0.485 0.6276
## x5 0.15467 0.02784 5.556 2.75e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.3264) family taken to be 1)
##
## Null deviance: 69.001 on 43 degrees of freedom
## Residual deviance: 42.648 on 38 degrees of freedom
## AIC: 319.85
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.3264
## Std. Err.: 0.0857
##
## 2 x log-likelihood: -305.8500#Comprobamos la calidad del ajuste basado en la deviance.
ajuste2 = 1-pchisq(42.648/0.3264,38)
ajuste2## [1] 4.262146e-12#Hay evidencia significativa que el ajuste del modelo no es apropiado (p<0.0001)
#Modelo sin intercepto
#fit3 <- glm.nb(y~-1+x1+x2+x3+x4+x5, data=datos)
#summary(fit3)
#fit4 <- glm.nb(y~x1+x2+x3+x5, data=datos)
#summary(fit4)
#Error en la estimación de este modelo, el exceso de cerosplot(table(datos$y))
Se observa que, efectivamente, una gran parte de los árboles no tienen plaga. En realidad, si los datos siguieran una distribución de Poisson, el número de ceros que cabría esperar es mucho menor que el observado en estos datos.
Para realizar los análisis exploratorios siguientes, utilizaremos el paquete pscl e inicialmente necesitaremos crear un objeto del tipo zeroinfl, como sigue.
library(pscl)## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim ZeileisModelos inflados por ceros
En ecología es comum encontrar en datos de conteos valores 0. Sin embargo, sí que es mucho más usual encontrar datos de conteos en los que hay un numero de ceros mayor que el que cabría esperar de acuerdo a una distribución de Poisson o una binomial negativa.
- Esto puede causar problemas en nuestros modelos ya que, de no tener en cuenta el exceso de ceros:
- Las estimas de los coeficientes pueden ser poco confiables.
- Puede haber sobredispersión.
Hay varios motivos por los quales se da la presecia ceros. Por ejemplo, en el contexto de espécies de pragas en un área forestal:
Hay errores estructurales. Es decir, que una plaga no está presente en un parche porque el hábitat no es adecuado.
Hay errores de diseño, debidos a un diseño experimental o muestral incorrecto. Por ejemplo, si buscamos a una especie en una época en la que los individuos se encuentran en otro sitio (e.g. espécies de plagas que atacam viveros florestais), es muy probable que nuestros conteos contengan una gran proporcion de ceros.
Hay errores de observador. Esto ocurre cuando dos especies son similares y el observador no sabe distinguirlas, o cuando son dificiles de detectar.
Nota
Asumir que los ceros proceden de dos procesos distintos: el proceso binomial y el proceso de Poisson. Igual que en los ZAP, se hace un glm binomial para modelar la probabilidad de medir un 0 (los falsos ceros; por ejemplo cuando el observador no ha visto ningún individuo de la especie en cuestión a pesar de que esté presente). Posteriormente se modela la probabilidad de obtener el resto de valores, incluyendo ceros. Estos modelos se llaman Modelos de mescla de distribuciones (mixture models), o modelos ZIP (Zero Inflated) o ZINB.
Definición...:
- Caso Poisson inflada por ceros: \(Y\sim ZIP(\mu_i,\pi_i)\)
Supongamos que \(y_i,\ldots,y_n\) son Observaciones de la variable respuesta \(Y_i,i=1,2,\ldots,n\), el modelo de Poisson inflado ceros se expresa como sigue,
\[ P(Y_i=y_i|x_i) =\left\{ \begin{array}{cccccccccccccccc} \{\pi_i+(1-\pi_i)e^{-\mu_i}, \;\; y_i=0 \\ (1-\pi_i)\frac{e^{-\mu_i}\mu_i^{y_i}} {y_i!}, \;\; y_i=1,2,3,\ldots,\\ \end{array} \right. \] por lo tanto \(0<\pi_i<1\) e \(\mu_i>0\).
La media y la varianza de la distribución son: \[ E[Y_i] = (1-\pi_i)\mu_i\\ Var[Y_i] = \mu_i(1-\pi_i)(1+\pi_i\mu_i) \]
El modelo ZIP, implica que la media \(\mu_i\) de una variable de Poisson mediante un modelo de regresión de Poisson y la probabilidad de \(\pi_i\) a través de una función de enlace de regresión logística \(\eta_i = logit(\pi_i)\).
- Caso Binomial Negativa inflada por ceros: \(Y\sim ZINB(\mu_i,\alpha,\pi_i)\)
El modelo Binomial Negativo inflado por ceros é expresado asi,
\[
P(Y_i=y_i|x_i) = \left\{
\begin{array}{ccccccc}
\pi_i +(1-\pi_i)\Biggl(\frac{1}{1+\alpha\mu}\Biggr)^{\alpha^{-1}} \;\; y_i=0 \\
(1-\pi_i)\frac{\Gamma(y_i+\alpha^{-1})}{y_i!\Gamma(\alpha^{-1})}\Biggl(\frac{\alpha\mu_i}{1+\alpha\mu_i}\Biggr)^{y_i}\Biggl(\frac{1}{1+\alpha\mu_i}\Biggr)^{\alpha^{-1}} \;\; y_i=1,2,\ldots\\
\end{array}
\right.
\] por lo tanto \(0<\pi_i<1\) , \(\mu_i>0\).
La media y la varianza de la distribución son:
\[
E[Y_i] = (1-\pi_i)\mu_i\\
Var[Y_i] = \mu_i(1-\pi_i)(1+\pi_i\mu_i+\alpha\mu_i)
\]
fit3 <- zeroinfl(y~x1+x2+x3+x4+x5, dist ="poisson",data=datos)
summary(fit3) ##
## Call:
## zeroinfl(formula = y ~ x1 + x2 + x3 + x4 + x5, data = datos, dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -5.6762 -0.5719 -0.1988 0.7403 7.9387
##
## Count model coefficients (poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.178669 1.549235 -0.115 0.90819
## x1 0.327276 0.048137 6.799 1.05e-11 ***
## x2 -0.245856 0.023178 -10.607 < 2e-16 ***
## x3 -0.235425 0.076722 -3.069 0.00215 **
## x4 0.857696 0.374711 2.289 0.02208 *
## x5 0.026471 0.003635 7.282 3.30e-13 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 51.68555 33.70215 1.534 0.12513
## x1 -1.44380 0.96677 -1.493 0.13533
## x2 -0.15271 0.48070 -0.318 0.75073
## x3 2.60999 1.62528 1.606 0.10830
## x4 -3.78812 6.56707 -0.577 0.56405
## x5 -0.19061 0.06525 -2.921 0.00349 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 1
## Log-likelihood: -288.8 on 12 Dffit4 <- zeroinfl(y~-1+x1+x2+x3+x4+x5, dist ="poisson",data=datos)
summary(fit4) ##
## Call:
## zeroinfl(formula = y ~ -1 + x1 + x2 + x3 + x4 + x5, data = datos, dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -5.6336 -0.7454 -0.2411 0.8147 6.9015
##
## Count model coefficients (poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## x1 0.322309 0.021487 15.00 < 2e-16 ***
## x2 -0.245304 0.022667 -10.82 < 2e-16 ***
## x3 -0.226755 0.015162 -14.96 < 2e-16 ***
## x4 0.816345 0.108551 7.52 5.46e-14 ***
## x5 0.026343 0.003462 7.61 2.74e-14 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## x1 -0.15391 0.38059 -0.404 0.68591
## x2 -0.06198 0.43523 -0.142 0.88676
## x3 0.21183 0.24276 0.873 0.38289
## x4 6.39928 2.46460 2.596 0.00942 **
## x5 -0.21895 0.07101 -3.083 0.00205 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 1
## Log-likelihood: -290.1 on 10 Dffit5 <- zeroinfl(y~-1+x1+x3+x4+x5, dist ="poisson",data=datos)
summary(fit5) ##
## Call:
## zeroinfl(formula = y ~ -1 + x1 + x3 + x4 + x5, data = datos, dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -7.5627 -0.6803 -0.2761 0.6929 4.6828
##
## Count model coefficients (poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## x1 0.093219 0.003568 26.125 < 2e-16 ***
## x3 -0.218265 0.015041 -14.511 < 2e-16 ***
## x4 1.077164 0.094980 11.341 < 2e-16 ***
## x5 0.021320 0.003268 6.524 6.85e-11 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## x1 -0.20657 0.09477 -2.180 0.02928 *
## x3 0.21338 0.24155 0.883 0.37703
## x4 6.27466 2.27483 2.758 0.00581 **
## x5 -0.21650 0.06826 -3.171 0.00152 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 1
## Log-likelihood: -353 on 8 Dffit6 <- zeroinfl(y~-1+x1+x4+x5, dist ="poisson",data=datos)
summary(fit6)##
## Call:
## zeroinfl(formula = y ~ -1 + x1 + x4 + x5, data = datos, dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -6.8455 -0.7447 -0.2769 0.6392 7.9356
##
## Count model coefficients (poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## x1 0.056449 0.002588 21.810 < 2e-16 ***
## x4 0.755010 0.093820 8.047 8.46e-16 ***
## x5 0.040632 0.003184 12.762 < 2e-16 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## x1 -0.15524 0.06753 -2.299 0.02152 *
## x4 5.91620 2.11126 2.802 0.00508 **
## x5 -0.21613 0.06792 -3.182 0.00146 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 1
## Log-likelihood: -461.4 on 6 DfNada de lo que hemos visto hasta ahora nos indica que nuestro modelo inflado por ceros sea mejor que una Poisson y corriente. El problema aquí es que es dificil comparar modelos que no están anidados. Sin embargo, el test de Vuong, disponible en la función pscl::vuong, nos permite hacer estas comparaciones.
vuong(fit6,fit1)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw 1.562327 model1 > model2 0.059106
## AIC-corrected 1.562327 model1 > model2 0.059106
## BIC-corrected 1.562327 model1 > model2 0.059106El test nos indica que el modelo inflado por ceros es mejor que su equivalente sin inflar.
Solo nos queda validar el modelo. Como siempre, calculamos los residuos de Pearson:
hnp(fit6, print=TRUE,verb=TRUE,xlab="Cuantiles teóricos") 
El modelo binomial negativo inflado de ceros - ZINB
fit7 <- zeroinfl(y~x1+x2+x3+x4+x5, dist ="negbin",data=datos)
summary(fit7) ##
## Call:
## zeroinfl(formula = y ~ x1 + x2 + x3 + x4 + x5, data = datos, dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -1.4080 -0.4350 -0.1675 0.3181 2.2677
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.92492 7.90167 0.117 0.90682
## x1 0.29284 0.21469 1.364 0.17256
## x2 -0.24109 0.11493 -2.098 0.03593 *
## x3 -0.14753 0.39225 -0.376 0.70684
## x4 0.45452 2.02726 0.224 0.82260
## x5 0.02927 0.01946 1.504 0.13246
## Log(theta) 0.86063 0.30349 2.836 0.00457 **
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 51.68560 33.73151 1.532 0.12546
## x1 -1.44277 0.96736 -1.491 0.13584
## x2 -0.15379 0.48128 -0.320 0.74932
## x3 2.60949 1.62658 1.604 0.10865
## x4 -3.78760 6.57424 -0.576 0.56453
## x5 -0.19075 0.06544 -2.915 0.00356 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 2.3646
## Number of iterations in BFGS optimization: 14
## Log-likelihood: -130.5 on 13 Dffit8 <- zeroinfl(y~-1+x1+x2+x3+x4+x5, dist ="negbin",data=datos)
summary(fit8)##
## Call:
## zeroinfl(formula = y ~ -1 + x1 + x2 + x3 + x4 + x5, data = datos, dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -1.3840 -0.4684 -0.2028 0.3628 1.8384
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## x1 0.31452 0.10813 2.909 0.00363 **
## x2 -0.23983 0.11446 -2.095 0.03614 *
## x3 -0.19241 0.08393 -2.293 0.02187 *
## x4 0.68243 0.56151 1.215 0.22423
## x5 0.02968 0.01918 1.548 0.12170
## Log(theta) 0.85803 0.30513 2.812 0.00492 **
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## x1 -0.15303 0.38110 -0.402 0.68802
## x2 -0.06308 0.43581 -0.145 0.88492
## x3 0.21125 0.24324 0.869 0.38512
## x4 6.40513 2.47062 2.593 0.00953 **
## x5 -0.21927 0.07143 -3.070 0.00214 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 2.3585
## Number of iterations in BFGS optimization: 16
## Log-likelihood: -131.8 on 11 Dffit8 <- zeroinfl(y~-1+x1+x3+x4+x5, dist ="negbin",data=datos)
summary(fit8)##
## Call:
## zeroinfl(formula = y ~ -1 + x1 + x3 + x4 + x5, data = datos, dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -1.2822 -0.4423 -0.1964 0.2802 1.7969
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## x1 0.09186 0.02203 4.170 3.05e-05 ***
## x3 -0.20910 0.09208 -2.271 0.0232 *
## x4 1.04021 0.58249 1.786 0.0741 .
## x5 0.02367 0.02109 1.122 0.2618
## Log(theta) 0.68494 0.30004 2.283 0.0224 *
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## x1 -0.20685 0.09512 -2.175 0.02966 *
## x3 0.21119 0.24268 0.870 0.38417
## x4 6.29650 2.29153 2.748 0.00600 **
## x5 -0.21766 0.06945 -3.134 0.00172 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 1.9836
## Number of iterations in BFGS optimization: 15
## Log-likelihood: -133.8 on 9 Dffit9 <- zeroinfl(y~-1+x1+x4+x5, dist ="negbin",data=datos)
summary(fit9)##
## Call:
## zeroinfl(formula = y ~ -1 + x1 + x4 + x5, data = datos, dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -1.2164 -0.4720 -0.2166 0.2068 1.8041
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## x1 0.05868 0.01849 3.173 0.00151 **
## x4 0.70406 0.65995 1.067 0.28605
## x5 0.04604 0.01968 2.339 0.01934 *
## Log(theta) 0.49148 0.29526 1.665 0.09600 .
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## x1 -0.15596 0.06813 -2.289 0.02207 *
## x4 5.93533 2.13370 2.782 0.00541 **
## x5 -0.21724 0.06955 -3.124 0.00179 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 1.6347
## Number of iterations in BFGS optimization: 12
## Log-likelihood: -136.5 on 7 DfEl test de Vuong, disponible en la función pscl::vuong, nos permite hacer estas comparaciones. Probemos a calcular el modelo binomial negativa (usando la función glm.nb de MASS):
vuong(fit9,fit2)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw 3.591767 model1 > model2 0.00016422
## AIC-corrected 3.591767 model1 > model2 0.00016422
## BIC-corrected 3.591767 model1 > model2 0.00016422Perfecto; el test nos indica que el modelo inflado por ceros es mucho mejor que su equivalente sin inflar.
Calculamos los residuos de Pearson:
hnp(fit9, print=TRUE,verb=TRUE,halfnormal =T,xlab="Cuantiles teóricos") 