quiz
GINNA VALENTINA ROMERO RINCON
2/11/2022
Para el conjunto de datos de Excel llamado poissonoffset.xlsx ajuste un modelo de Poisson para las variables edad de la planta(años)[ manejarla de modo numérica (no factor)] y para el tipo de híbrido usando como respuesta el número de abortos contados para ese año. Existen dos variables llamadas plantas y plantas2 que nos sirven para ajustar el modelo en dos escenarios, uno para ajustar offset pues el número de plantas contadas fue diferente y el otro fue el mismo en todos los muestreos. El offset nos permite ajustar tasas y no solo conteos.
Importando los datos
#Se importan los datos y se cambia a numeric la variable abortos
library(readxl)
datos1 <- read_excel("C:/Users/ginna/Desktop/ModelacionCultivos/quiz/datos1.xlsx", col_types = c("numeric", "numeric", "numeric",
"numeric", "text", "numeric", "numeric"))
## creando el data frame
dfdatos = data.frame(datos1)• Involucre en el modelo la edad y la edad al cuadrado • Ajuste ambos modelos (con y sin offset) • Escriba ambas ecuaciones
#Para el modelo poisson 1 que se crea se utiliza offset para estandarizar los datos ya que el numero de plantas muestreado es desigual en los conteos
mod1<- glm(abortos ~ edad + I(edad^2)
+ hibrido + (hibrido*edad) + (hibrido*I(edad^2))
+ offset(log(plantas)), family = poisson(link = "log"),
data=dfdatos)
summary(mod1)##
## Call:
## glm(formula = abortos ~ edad + I(edad^2) + hibrido + (hibrido *
## edad) + (hibrido * I(edad^2)) + offset(log(plantas)), family = poisson(link = "log"),
## data = dfdatos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.59230 -0.29255 -0.04407 0.28813 0.67432
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.660908 1.461440 -1.821 0.0686 .
## edad 0.258879 0.269497 0.961 0.3368
## I(edad^2) -0.002283 0.011917 -0.192 0.8481
## hibridoh2 -0.189849 2.165070 -0.088 0.9301
## edad:hibridoh2 0.050225 0.396211 0.127 0.8991
## I(edad^2):hibridoh2 -0.001993 0.017419 -0.114 0.9089
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 87.5202 on 47 degrees of freedom
## Residual deviance: 6.2911 on 42 degrees of freedom
## AIC: 171.41
##
## Number of Fisher Scoring iterations: 4Ecuacion Modelo 1
\[log() = -2,7 + 0,3(edad) -0,002(edad^{2}) - 0,3(h2) + 0,07(edad:h2) - 0,003 (edad^{2}:h2)\]
#Para el modelo poisson 2 se elimina offset
mod2<- glm(abortos ~ edad + I(edad^2)
+ hibrido + (hibrido*edad) + (hibrido*I(edad^2)),
family = poisson(link = "log"),
data=dfdatos)
summary(mod2)##
## Call:
## glm(formula = abortos ~ edad + I(edad^2) + hibrido + (hibrido *
## edad) + (hibrido * I(edad^2)), family = poisson(link = "log"),
## data = dfdatos)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.44194 -0.26641 0.08525 0.27340 0.46031
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.936936 1.453066 -0.645 0.519
## edad 0.229571 0.268007 0.857 0.392
## I(edad^2) -0.001067 0.011856 -0.090 0.928
## hibridoh2 -0.628668 2.159554 -0.291 0.771
## edad:hibridoh2 0.131838 0.395308 0.334 0.739
## I(edad^2):hibridoh2 -0.005984 0.017385 -0.344 0.731
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 80.2678 on 47 degrees of freedom
## Residual deviance: 3.7591 on 42 degrees of freedom
## AIC: 168.88
##
## Number of Fisher Scoring iterations: 4Ecuacion Modelo 2
\[log() = -0,9 + 0,2(edad) -0,001(edad^{2}) - 0,7(h2) + 0,1(edad:h2) - 0,006 (edad^{2}:h2)\]
# Para realizar la comparacion entre los modelos se utiliza la funcion cbind
library(arm)## Loading required package: MASS## Loading required package: Matrix## Loading required package: lme4## Warning: package 'lme4' was built under R version 4.0.5##
## arm (Version 1.13-1, built: 2022-8-25)## Working directory is C:/Users/ginna/Desktopcoef1 = coef(mod1)
# extract coefficients from second model
coef2 = coef(mod2)
# extract standard errors from first model using 'se.coef()'
se.coef1 = se.coef(mod1)
# extract standard errors from second model
se.coef2 = se.coef(mod2)
# use 'cbind()' to combine values into one dataframe
models.both <- cbind(coef1, se.coef1, coef2, se.coef2, exponent = exp(coef1))
# show dataframe
models.both## coef1 se.coef1 coef2 se.coef2 exponent
## (Intercept) -2.660908161 1.46144025 -0.936936138 1.45306584 0.06988473
## edad 0.258879044 0.26949727 0.229570540 0.26800666 1.29547710
## I(edad^2) -0.002282629 0.01191696 -0.001067061 0.01185573 0.99771997
## hibridoh2 -0.189848919 2.16507004 -0.628667830 2.15955449 0.82708408
## edad:hibridoh2 0.050224824 0.39621063 0.131838095 0.39530774 1.05150747
## I(edad^2):hibridoh2 -0.001993064 0.01741864 -0.005984131 0.01738457 0.99800892• Juzgue la dispersión del modelo • Seleccione el modelo que no muestre problemas de sobre o sub-dispersión
En la prueba de dispersiontest se evalua si nos da un valor >1 puede ser sobredispersion; si es < 1 entonces es subdispersion y si es igual entonces no presenta dispercion
library(AER)## Loading required package: car## Warning: package 'car' was built under R version 4.0.5## Loading required package: carData## Warning: package 'carData' was built under R version 4.0.3##
## Attaching package: 'car'## The following object is masked from 'package:arm':
##
## logit## Loading required package: lmtest## Warning: package 'lmtest' was built under R version 4.0.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.0.5##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Loading required package: survivaldt1 <- dispersiontest(mod1,alternative = "t")
dt1##
## Dispersion test
##
## data: mod1
## z = -27.935, p-value < 2.2e-16
## alternative hypothesis: true dispersion is not equal to 1
## sample estimates:
## dispersion
## 0.1194403dt2 <- dispersiontest(mod2,alternative = "t")
dt2##
## Dispersion test
##
## data: mod2
## z = -35.003, p-value < 2.2e-16
## alternative hypothesis: true dispersion is not equal to 1
## sample estimates:
## dispersion
## 0.07902589Para ambos modelos la dispersion no es la mas adecuada sin embargo entre los dos modelos se selecciona el modelo 1 donde se utiliza offset
• Elimine términos “no significativos” para reducir el modelo y demuestre si este modelo reducido resulta más conveniente que el modelo ampliado (función step para quitar mugre)
stats::step(mod1)## Start: AIC=171.41
## abortos ~ edad + I(edad^2) + hibrido + (hibrido * edad) + (hibrido *
## I(edad^2)) + offset(log(plantas))
##
## Df Deviance AIC
## - I(edad^2):hibrido 1 6.3042 169.43
## - edad:hibrido 1 6.3071 169.43
## <none> 6.2911 171.41
##
## Step: AIC=169.43
## abortos ~ edad + I(edad^2) + hibrido + edad:hibrido + offset(log(plantas))
##
## Df Deviance AIC
## - edad:hibrido 1 6.3154 167.44
## - I(edad^2) 1 6.4421 167.56
## <none> 6.3042 169.43
##
## Step: AIC=167.44
## abortos ~ edad + I(edad^2) + hibrido + offset(log(plantas))
##
## Df Deviance AIC
## - I(edad^2) 1 6.4508 165.57
## - hibrido 1 6.9761 166.10
## <none> 6.3154 167.44
## - edad 1 8.4160 167.54
##
## Step: AIC=165.57
## abortos ~ edad + hibrido + offset(log(plantas))
##
## Df Deviance AIC
## - hibrido 1 7.110 164.23
## <none> 6.451 165.57
## - edad 1 86.293 243.41
##
## Step: AIC=164.23
## abortos ~ edad + offset(log(plantas))
##
## Df Deviance AIC
## <none> 7.11 164.23
## - edad 1 87.52 242.64##
## Call: glm(formula = abortos ~ edad + offset(log(plantas)), family = poisson(link = "log"),
## data = dfdatos)
##
## Coefficients:
## (Intercept) edad
## -2.3747 0.2105
##
## Degrees of Freedom: 47 Total (i.e. Null); 46 Residual
## Null Deviance: 87.52
## Residual Deviance: 7.11 AIC: 164.2Al evaluar las variables que deberian quedarse se selecciona las siguientes variables
\[glm(formula abortos ~ edad + offset(log(plantas)), family = poisson(link = "log"), data = dfdatos)\]
por lo tanto el modelo quedaria de la siguiente manera:
modfinal<- glm(abortos ~ edad,
offset(log(plantas)),
family = poisson(),
data=dfdatos)
summary(modfinal)##
## Call:
## glm(formula = abortos ~ edad, family = poisson(), data = dfdatos,
## weights = offset(log(plantas)))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.61440 -0.34189 0.04171 0.35003 0.50880
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.78063 0.24283 -3.215 0.00131 **
## edad 0.20480 0.01973 10.381 < 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: 122.5647 on 47 degrees of freedom
## Residual deviance: 6.4536 on 46 degrees of freedom
## AIC: 244.12
##
## Number of Fisher Scoring iterations: 4\[log() = -2,3 + 0,2(edad) \]
• Estime los conteos esperados para las diferentes edades del híbrido 2 (dar valores)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.0.5dfdatos$pred=fitted.values(modfinal)
g1 <- ggplot(dfdatos, aes(abortos, pred, color = hibrido))+
geom_point()+
theme_classic()
g1
g1 <- ggplot(dfdatos, aes(edad, pred, color = hibrido))+
geom_point(size = (dfdatos$edad))+
theme_classic()
g1
¿percibe alguna tendencia que sugiera alguna relación clara con la edad? Si parece estar relacionado que a mayor edad mayor es el numero de abortos