Quiz modelación de cultivos
Luisa Marina Reyes Forero 1022407355
2022-11-02
#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.
library(readxl)
poisson <- read_excel("C:/Users/AGUA MARINA/Downloads/poisson.xlsx")
View(poisson)#Involucre en el modelo la edad y la edad al cuadrado #Ajuste ambos modelos (con y sin offset) #Escriba ambas ecuaciones
#modelo sin offset
abortossin=glm(abortos~edad+edad2+hibrido+edad*hibrido+edad2*hibrido, family = poisson(),data = poisson)
summary(abortossin)##
## Call:
## glm(formula = abortos ~ edad + edad2 + hibrido + edad * hibrido +
## edad2 * hibrido, family = poisson(), data = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5904 -0.2728 -0.1237 0.3058 1.1777
##
## 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
## edad2 -0.001067 0.011856 -0.090 0.928
## hibridoh2 -0.299276 2.146618 -0.139 0.889
## edad:hibridoh2 0.055696 0.392596 0.142 0.887
## edad2:hibridoh2 -0.001844 0.017241 -0.107 0.915
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 90.3020 on 47 degrees of freedom
## Residual deviance: 5.6536 on 42 degrees of freedom
## AIC: 171.18
##
## Number of Fisher Scoring iterations: 4library(sjPlot)## Warning: package 'sjPlot' was built under R version 4.2.2## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!tab_model(abortossin,
transform = NULL,
show.r2 = FALSE)| abortos | |||
|---|---|---|---|
| Predictors | Log-Mean | CI | p |
| (Intercept) | -0.94 | -3.93 – 1.79 | 0.519 |
| edad | 0.23 | -0.28 – 0.77 | 0.392 |
| edad2 | -0.00 | -0.02 – 0.02 | 0.928 |
| hibrido [h2] | -0.30 | -4.56 – 3.90 | 0.889 |
| edad * hibrido [h2] | 0.06 | -0.71 – 0.83 | 0.887 |
| edad2 * hibrido [h2] | -0.00 | -0.04 – 0.03 | 0.915 |
| Observations | 48 | ||
#modelo con offset
abortoscon=glm(abortos~edad+edad2+hibrido+edad*hibrido+edad2*hibrido, family = poisson(),data = poisson,offset = plantas)
summary(abortoscon)##
## Call:
## glm(formula = abortos ~ edad + edad2 + hibrido + edad * hibrido +
## edad2 * hibrido, family = poisson(), data = poisson, offset = plantas)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6383 -0.5843 -0.1687 1.0340 2.3683
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.649550 1.527820 -4.352 1.35e-05 ***
## edad 0.369744 0.280937 1.316 0.188
## edad2 -0.006441 0.012373 -0.521 0.603
## hibridoh2 1.743253 2.195735 0.794 0.427
## edad:hibridoh2 -0.311417 0.400888 -0.777 0.437
## edad2:hibridoh2 0.015627 0.017564 0.890 0.374
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 174.983 on 47 degrees of freedom
## Residual deviance: 56.238 on 42 degrees of freedom
## AIC: 221.76
##
## Number of Fisher Scoring iterations: 5library(sjPlot)
tab_model(abortoscon,
transform = NULL,
show.r2 = FALSE)| abortos | |||
|---|---|---|---|
| Predictors | Log-Mean | CI | p |
| (Intercept) | -6.65 | -9.79 – -3.78 | <0.001 |
| edad | 0.37 | -0.17 – 0.94 | 0.188 |
| edad2 | -0.01 | -0.03 – 0.02 | 0.603 |
| hibrido [h2] | 1.74 | -2.61 – 6.04 | 0.427 |
| edad * hibrido [h2] | -0.31 | -1.10 – 0.48 | 0.437 |
| edad2 * hibrido [h2] | 0.02 | -0.02 – 0.05 | 0.374 |
| Observations | 48 | ||
#ecuación de modelo abortos sin offset \[abortossin=0.39+1.26edad+1edad^2+0.74hibrido+1.06edad*hibrido+1edad^2*hibrido\] #ecuació de modelo con offset \[abortoscon=-6.65+0.37edad-0.01edad^2+1.74hibrido-0.31edad*hibrido+0.02edad^2*hibrido\] #Juzgue la dispersión del modelo #El modelo sin ajustar presenta una subdispersión lo que lleva a pensar en que su varianza es mucho menor respecto a la esperada, adicionalmente, los errores estándar estimados y los estadísticos de prueba de la bondad de ajuste general están distorsionados, por tanto si hay que realizar un ajuste. #El modelo ajustado no presenta dispersión, por tanto el ajuste que se le realizo fue el adecuado.
#Seleccione el modelo que no muestre problemas de sobre o sub-dispersión
library(AER)## Warning: package 'AER' was built under R version 4.2.2## Loading required package: car## Warning: package 'car' was built under R version 4.2.2## Loading required package: carData## Warning: package 'carData' was built under R version 4.2.2## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Warning: package 'sandwich' was built under R version 4.2.2## Loading required package: survivaldispersiontest(abortossin,alternative = "t")##
## Dispersion test
##
## data: abortossin
## z = -24.56, p-value < 2.2e-16
## alternative hypothesis: true dispersion is not equal to 1
## sample estimates:
## dispersion
## 0.1212782dispersiontest(abortoscon,alternative = "t")##
## Dispersion test
##
## data: abortoscon
## z = 0.34922, p-value = 0.7269
## alternative hypothesis: true dispersion is not equal to 1
## sample estimates:
## dispersion
## 1.080064#Elimine términos “no significativos” para reducir el modelo y demuestre si este modelo reducido resulta más conveniente que el modelo ampliado
step(abortoscon)## Start: AIC=221.76
## abortos ~ edad + edad2 + hibrido + edad * hibrido + edad2 * hibrido
##
## Df Deviance AIC
## - edad:hibrido 1 56.840 220.36
## - edad2:hibrido 1 57.028 220.55
## <none> 56.238 221.76
##
## Step: AIC=220.36
## abortos ~ edad + edad2 + hibrido + edad2:hibrido
##
## Df Deviance AIC
## - edad2:hibrido 1 57.789 219.31
## - edad 1 58.075 219.60
## <none> 56.840 220.36
##
## Step: AIC=219.31
## abortos ~ edad + edad2 + hibrido
##
## Df Deviance AIC
## - edad2 1 57.821 217.34
## - edad 1 58.906 218.43
## <none> 57.789 219.31
## - hibrido 1 65.715 225.24
##
## Step: AIC=217.34
## abortos ~ edad + hibrido
##
## Df Deviance AIC
## <none> 57.821 217.34
## - hibrido 1 65.717 223.24
## - edad 1 165.784 323.31##
## Call: glm(formula = abortos ~ edad + hibrido, family = poisson(), data = poisson,
## offset = plantas)
##
## Coefficients:
## (Intercept) edad hibridoh2
## -6.1273 0.2447 0.3795
##
## Degrees of Freedom: 47 Total (i.e. Null); 45 Residual
## Null Deviance: 175
## Residual Deviance: 57.82 AIC: 217.3#el mejor modelo es el de abortos con ajuste ya que no presento problemas de dispersión al obtener 1 en la prueba de dispersión de Poisson. Se toma la decision de no trabajar con modelo quasi-poisson por esá razón.
abortosconmejor=glm(abortos~edad+hibrido, family = poisson(),data = poisson,offset = plantas)
summary(abortoscon)##
## Call:
## glm(formula = abortos ~ edad + edad2 + hibrido + edad * hibrido +
## edad2 * hibrido, family = poisson(), data = poisson, offset = plantas)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6383 -0.5843 -0.1687 1.0340 2.3683
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.649550 1.527820 -4.352 1.35e-05 ***
## edad 0.369744 0.280937 1.316 0.188
## edad2 -0.006441 0.012373 -0.521 0.603
## hibridoh2 1.743253 2.195735 0.794 0.427
## edad:hibridoh2 -0.311417 0.400888 -0.777 0.437
## edad2:hibridoh2 0.015627 0.017564 0.890 0.374
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 174.983 on 47 degrees of freedom
## Residual deviance: 56.238 on 42 degrees of freedom
## AIC: 221.76
##
## Number of Fisher Scoring iterations: 5library(sjPlot)
tab_model(abortosconmejor,
transform = NULL,
show.r2 = FALSE)| abortos | |||
|---|---|---|---|
| Predictors | Log-Mean | CI | p |
| (Intercept) | -6.13 | -6.75 – -5.53 | <0.001 |
| edad | 0.24 | 0.20 – 0.29 | <0.001 |
| hibrido [h2] | 0.38 | 0.12 – 0.64 | 0.005 |
| Observations | 48 | ||
#ecuación modelo reducido \[abortosconreducido=-6.13+0.24edad+0.38hibrido\]
library(AER)
dispersiontest(abortosconmejor,alternative = "t")##
## Dispersion test
##
## data: abortosconmejor
## z = 0.52088, p-value = 0.6024
## alternative hypothesis: true dispersion is not equal to 1
## sample estimates:
## dispersion
## 1.120911#Estime los conteos esperados para las diferentes edades del híbrido 2 #Grafique los conteos esperados como función de la edad de la planta en el híbrido 2. ¿percibe alguna tendencia que sugiera alguna relación clara con la edad? #si se percibe una realción entre los conteos esperados y la edad ya que a medida que va teniendo mas años el conteo de abortos aumenta para el hibrido 2 por tanto esta es una señal de que la edad si influye en este comportamiento de la planta. Es muy poca la diferencia ya que es en decimales por tanto este dato tambien puede ser despreciable por la cantidad de años que tieen que pasar para que sea una doferencia realemnte marcada.
library(sjPlot)
hibrido2=plot_model(abortoscon,"pred",
terms=c("edad","hibrido"),
axis.labels = c("plantas","abortos"),
title = "Predicción abortos")
hibrido2$data## # Predicted counts of abortos
##
## # hibrido = h1
##
## edad | Predicted | group_col | 95% CI
## ---------------------------------------------
## 6 | 0.59 | h1 | [0.05, 6.60]
## 8 | 1.25 | h1 | [0.34, 4.63]
## 10 | 2.61 | h1 | [1.96, 3.48]
## 12 | 5.47 | h1 | [2.15, 13.88]
## 15 | 16.58 | h1 | [1.27, 216.07]
##
## # hibrido = h2
##
## edad | Predicted | group_col | 95% CI
## --------------------------------------------
## 6 | 3.12 | h2 | [0.27, 36.48]
## 8 | 3.50 | h2 | [0.91, 13.48]
## 10 | 3.94 | h2 | [2.89, 5.36]
## 12 | 4.42 | h2 | [1.72, 11.34]
## 15 | 5.27 | h2 | [0.39, 71.29]
##
## Adjusted for:
## * edad2 = 114.00
## * plantas = 4.65hibrido2
#¿cuantas veces es mayor el conteo de un híbrido respecto al otro para la edad de 15
hibrido2=data.frame(poisson$abortos,poisson$edad,poisson$hibrido)
hibrido2## poisson.abortos poisson.edad poisson.hibrido
## 1 2 6 h1
## 2 2 6 h1
## 3 1 6 h1
## 4 2 6 h1
## 5 1 6 h1
## 6 1 6 h2
## 7 1 6 h2
## 8 2 6 h2
## 9 1 6 h2
## 10 3 8 h1
## 11 2 8 h1
## 12 2 8 h1
## 13 2 8 h1
## 14 2 8 h1
## 15 2 8 h2
## 16 3 8 h2
## 17 3 8 h2
## 18 3 8 h2
## 19 3 10 h1
## 20 3 10 h1
## 21 4 10 h1
## 22 3 10 h1
## 23 3 10 h1
## 24 4 10 h2
## 25 4 10 h2
## 26 4 10 h2
## 27 3 10 h2
## 28 6 12 h1
## 29 6 12 h1
## 30 5 12 h1
## 31 5 12 h1
## 32 6 12 h1
## 33 6 12 h2
## 34 5 12 h2
## 35 6 12 h2
## 36 6 12 h2
## 37 6 12 h1
## 38 5 12 h2
## 39 5 12 h1
## 40 6 12 h2
## 41 9 15 h1
## 42 10 15 h1
## 43 10 15 h1
## 44 9 15 h2
## 45 10 15 h2
## 46 10 15 h2
## 47 9 15 h1
## 48 15 15 h2l2=hibrido2[41:48,]
l2## poisson.abortos poisson.edad poisson.hibrido
## 41 9 15 h1
## 42 10 15 h1
## 43 10 15 h1
## 44 9 15 h2
## 45 10 15 h2
## 46 10 15 h2
## 47 9 15 h1
## 48 15 15 h2l3H1=l2[c(1:3,7),]
l3H1## poisson.abortos poisson.edad poisson.hibrido
## 41 9 15 h1
## 42 10 15 h1
## 43 10 15 h1
## 47 9 15 h1library(dplyr)## Warning: package 'dplyr' was built under R version 4.2.2##
## 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, unionTotPesoshibrido1 = colSums (select (l3H1, contains ("abortos")))
TotPesoshibrido1## poisson.abortos
## 38l3H2=l2[c(4:6,8),]
l3H2## poisson.abortos poisson.edad poisson.hibrido
## 44 9 15 h2
## 45 10 15 h2
## 46 10 15 h2
## 48 15 15 h2library(dplyr)
TotPesoshibrido2 = colSums (select (l3H2, contains ("abortos")))
TotPesoshibrido2## poisson.abortos
## 44TASACONTEOS=TotPesoshibrido1/TotPesoshibrido2
TASACONTEOS## poisson.abortos
## 0.8636364#Es 0,86 veces mayor el conteo de abortos del hibrido 2 respecto al conteo de abortos del hibrido 1 para la edad 15.