Modelos de regresión de Poisson y binomial negativo
Data Fumadores
Fumadores.data <- read.table(file="http://tarwi.lamolina.edu.pe/~clopez/Categoricos/Fumadores.txt", header=T)
Fumadores.data[,1] <- c(39.5, 48.5, 57.5, 66.5, 75.5, 39.5, 48.5, 57.5, 66.5, 75.5)
attach(Fumadores.data)
Fumadores.data
## Edad Fumador Muertes Personas.Año
## 1 39.5 Si 32 52407
## 2 48.5 Si 104 43248
## 3 57.5 Si 206 28612
## 4 66.5 Si 186 12663
## 5 75.5 Si 102 5317
## 6 39.5 No 2 18790
## 7 48.5 No 12 10673
## 8 57.5 No 28 5710
## 9 66.5 No 28 2585
## 10 75.5 No 31 1462
plot(Edad, Muertes/Personas.Año)
m1 <- glm(Muertes ~ Edad + Fumador, family=poisson(link=log), offset=log(Personas.Año))
summary(m1)
##
## Call:
## glm(formula = Muertes ~ Edad + Fumador, family = poisson(link = log),
## offset = log(Personas.Año))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.5712 -2.7562 0.2857 1.4261 3.7183
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.950876 0.219729 -49.838 < 2e-16 ***
## Edad 0.092870 0.003227 28.777 < 2e-16 ***
## FumadorSi 0.406370 0.107195 3.791 0.00015 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 935.067 on 9 degrees of freedom
## Residual deviance: 69.182 on 7 degrees of freedom
## AIC: 130.25
##
## Number of Fisher Scoring iterations: 4
Residuales de Pearson y Devianza
res.pearson <- resid(m1, type="pearson")
res.est.p <- res.pearson/sqrt(1-lm.influence(m1)$hat)
res.devianza <- resid(m1, type="deviance")
res.est.d <- res.devianza/sqrt(1-lm.influence(m1)$hat)
residuales <- cbind(res.pearson,res.est.p,res.devianza,res.est.d)
residuales
## res.pearson res.est.p res.devianza res.est.d
## 1 -3.0036894 -3.7299251 -3.2536856 -4.0403657
## 2 0.1017365 0.1249275 0.1015672 0.1247196
## 3 3.8979086 4.5684480 3.7183048 4.3579477
## 4 2.0193510 2.4158510 1.9689713 2.3555793
## 5 -4.2823432 -6.5833379 -4.5712427 -7.0274693
## 6 -3.0376276 -3.3717960 -3.7912133 -4.2082834
## 7 -1.1971988 -1.3333823 -1.2637026 -1.4074511
## 8 1.5562726 1.7498813 1.4785350 1.6624727
## 9 1.3248189 1.5035311 1.2686020 1.4397308
## 10 0.4766739 0.5897146 0.4698269 0.5812438
Estadístico Chi-cuadrado y Devianza
X2 <- sum(res.pearson^2)
X2
## [1] 61.70725
p.valor <- 1-pchisq(q=X2, df=7)
p.valor
## [1] 6.879752e-11
D <- sum(res.devianza^2)
D
## [1] 69.18208
p.valor <- 1-pchisq(q=D, df=7)
p.valor
## [1] 2.161382e-12
Estadistico C y Pseudo R2
m0 <- glm(Muertes ~ 1, family=poisson(link=log))
C <- 2*(logLik(m1)[1] - logLik(m0)[1])
C
## [1] 575.0869
p.valor <- 1-pchisq(q=C, df=2)
p.valor
## [1] 0
pseudoR2 <- (logLik(m0)[1] - logLik(m1)[1])/logLik(m0)[1]
pseudoR2
## [1] 0.8223317
Regresión de Poisson
Data_Cangrejo <- read.table(file="http://tarwi.lamolina.edu.pe/~clopez/Categoricos/Cangrejo.txt", header=T)
head(Data_Cangrejo)
## Color Cond Ancho Sat Peso
## 1 3 3 28.3 8 3050
## 2 4 3 22.5 0 1550
## 3 2 1 26.0 9 2300
## 4 4 3 24.8 0 2100
## 5 4 3 26.0 4 2600
## 6 3 3 23.8 0 2100
attach(Data_Cangrejo)
Cangrejo.m1 <- glm(Sat ~ Peso*Ancho, family=poisson(link="log"))
summary(Cangrejo.m1)
##
## Call:
## glm(formula = Sat ~ Peso * Ancho, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1121 -1.8424 -0.5581 0.9171 4.9419
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.911e+00 1.906e+00 -3.626 0.000288 ***
## Peso 3.066e-03 7.467e-04 4.107 4.01e-05 ***
## Ancho 2.303e-01 7.520e-02 3.063 0.002191 **
## Peso:Ancho -8.661e-05 2.479e-05 -3.493 0.000477 ***
## ---
## 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: 545.29 on 169 degrees of freedom
## AIC: 908.59
##
## Number of Fisher Scoring iterations: 6
Test de Bohning
library(epiR)
## Loading required package: survival
## Loading required package: splines
## Package epiR 0.9-62 is loaded
## Type help(epi.about) for summary information
esperado <- predict(Cangrejo.m1, type="response")
epi.bohning(Sat, esperado)
## test.statistic p.value
## 1 19.59947 0
Regresión binomial negativa
library(MASS)
Cangrejo.NB1 <- glm(Sat ~ Peso*Ancho, family=negative.binomial(theta=1.0, link="log"), start=coef(Cangrejo.m1))
summary(Cangrejo.NB1)
##
## Call:
## glm(formula = Sat ~ Peso * Ancho, family = negative.binomial(theta = 1,
## link = "log"), start = coef(Cangrejo.m1))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8923 -1.3683 -0.2517 0.4127 2.2202
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.505e+00 3.676e+00 -2.042 0.0427 *
## Peso 3.416e-03 1.490e-03 2.293 0.0231 *
## Ancho 2.460e-01 1.469e-01 1.674 0.0959 .
## Peso:Ancho -9.695e-05 5.081e-05 -1.908 0.0581 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1) family taken to be 0.9308966)
##
## Null deviance: 224.93 on 172 degrees of freedom
## Residual deviance: 200.32 on 169 degrees of freedom
## AIC: 753.51
##
## Number of Fisher Scoring iterations: 4