Punto 1

Los siguientes datos corresponden a los pesos al nacer (en gramos) y la edad gestacional estimada (en semanas) de 24 bebés.

datos_nacimiento <- data.frame(
  Sexo = c(
    rep("Mujer", 12),
    rep("Hombre", 12)
  ),
  Edad = c(
    40, 36, 40, 38, 42, 39, 40, 37, 36, 38, 39, 40,
    40, 38, 40, 35, 36, 37, 41, 40, 37, 38, 40, 38
  ),
  Peso = c(
    3317, 2729, 2935, 2754, 3210, 2817, 3126, 2539, 2412, 2991, 2875, 3231,
    2968, 2795, 3163, 2925, 2625, 2847, 3292, 3473, 2628, 3176, 3421, 2975
  )
)
  1. Ajuste un modelo de respuesta normal con función de enlace identidad y las covariables edad gestacional, género y la interacción entre estas.

\[\begin{align*} \textbf{} & \quad Y_i \sim N(\mu_i, \phi) \\ \textbf{} & \quad g(\mu_i) = \eta_i = \mu_i \quad (\text{Identidad}) \\ \textbf{} & \quad \eta_i = \beta_0 + \beta_1 \text{Edad}_i + \beta_2 \text{Sexo}_i + \beta_3 (\text{Edad}_i \times \text{Sexo}_i) \\ \end{align*}\]

modelo_1.1 <- glm(
  Peso ~ Edad * Sexo,
  family = gaussian(link = "identity"),
  data = datos_nacimiento
)
  1. Evalúe, al nivel de significancia del 5%, si el efecto de la edad gestacional sobre el peso al nacer es el mismo para niños y niñas.

La covariable sexo tomará el valor de 0 cuando el \(i\)-ésimo bebé sea hombre y 1 cuando sea mujer.

  • Si es hombre (\(\text{Sexo}=0\)): \(\mu_i = \beta_0 + \beta_1 \text{Edad}_i\)
  • Si es mujer (\(\text{Sexo}=1\)): \(\mu_i = (\beta_0 + \beta_2) + (\beta_1+ \beta_3) \text{Edad}_i\)

Para que el efecto sea el mismo debe \(\beta_3 = 0\).

\[H_0: \beta_3 = 0 \quad \text{vs} \quad H_1: \beta_3 \neq 0\]

summary(modelo_1.1)
## 
## Call:
## glm(formula = Peso ~ Edad * Sexo, family = gaussian(link = "identity"), 
##     data = datos_nacimiento)
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -1268.67    1114.64  -1.138 0.268492    
## Edad             111.98      29.05   3.855 0.000986 ***
## SexoMujer       -872.99    1611.33  -0.542 0.593952    
## Edad:SexoMujer    18.42      41.76   0.441 0.663893    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 32621.23)
## 
##     Null deviance: 1829873  on 23  degrees of freedom
## Residual deviance:  652425  on 20  degrees of freedom
## AIC: 323.16
## 
## Number of Fisher Scoring iterations: 2

Interpretación: Ya que el valor-p del coeficiente de la interacción es mayor al nivel de significancia \(\alpha=0.05\) , no hay evidencia estadisticamente significativa para rechazar la hipótesis nula por lo que el efecto de la edad es el mismo en ambos sexos.

Se ajustara el modelo reducido (sin interacción) que evita problemas de multicolinealidad:

modelo_1.2 <- glm(
  Peso ~ Edad + Sexo,
  family = gaussian(link = "identity"),
  data = datos_nacimiento
)

summary(modelo_1.2)
## 
## Call:
## glm(formula = Peso ~ Edad + Sexo, family = gaussian(link = "identity"), 
##     data = datos_nacimiento)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1610.28     786.08  -2.049   0.0532 .  
## Edad          120.89      20.46   5.908 7.28e-06 ***
## SexoMujer    -163.04      72.81  -2.239   0.0361 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 31370.04)
## 
##     Null deviance: 1829873  on 23  degrees of freedom
## Residual deviance:  658771  on 21  degrees of freedom
## AIC: 321.39
## 
## Number of Fisher Scoring iterations: 2
  1. Realice el análisis de diagnóstico. Comente.
library(performance)
check_model(modelo_1.2)

Se observa que los valores predichos por el modelo se encuentran bastante próximos a los valores realmente observados, lo cual sugiere un buen ajuste general. Sin embargo, en lo referente a la linealidad, sí se evidencian problemas: la línea verde, que idealmente debería ser horizontal, muestra curvaturas que indican que la función de enlace podria no está logrando linealizar adecuadamente la relación entre las covariables y la media.

Por otro lado, no se detecta evidencia de multicolinealidad y sí aparecen algunos puntos que podrían considerarse influyentes, lo cual lleva a examinar con más detalle dichas observaciones.

Gráfico de envolvente (simulación):

El gráfico muestra un comportamiento satisfactorio, ya que los residuales se mantienen dentro de las bandas de confianza. Esto confirma que la distribución asumida por el modelo es apropiada. Dado que en el paso anterior también se descartó la presencia de multicolinealidad, se puede afirmar que el modelo ajusta adecuadamente a los datos.

Distancia de Cook (valores influyentes):

plot(modelo_1.2, which = 4, main = "Distancia de Cook")

Al analizar la influencia de las observaciones, se destaca el dato número 16, cuya distancia de Cook es cercana a 0.4. Para evaluar el impacto de este valor en el modelo, se comprobara si su eliminación modificaría de manera importante las estimaciones. En caso de que produzca cambios drásticos, es conveniente considerar un modelo que excluya esta observación.

modelo_1.2.1 <- glm(
  Peso ~ Edad + Sexo,
  family = gaussian(link = "identity"),
  data = datos_nacimiento[-16, ]
)

summary(modelo_1.2.1)
## 
## Call:
## glm(formula = Peso ~ Edad + Sexo, family = gaussian(link = "identity"), 
##     data = datos_nacimiento[-16, ])
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2318.03     801.57  -2.892  0.00902 ** 
## Edad          138.50      20.71   6.688 1.65e-06 ***
## SexoMujer    -137.40      68.54  -2.005  0.05870 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 26925.39)
## 
##     Null deviance: 1827974  on 22  degrees of freedom
## Residual deviance:  538508  on 20  degrees of freedom
## AIC: 304.68
## 
## Number of Fisher Scoring iterations: 2

Dado que no se observan cambios drásticos en las estimaciones de los coeficientes, se considera mas apropiado trabajar con el modelo que incluye todos los datos.

  1. Interprete las estimaciones de los parámetros, excluyendo el intercepto.
modelo_1.2$coefficients
## (Intercept)        Edad   SexoMujer 
##  -1610.2825    120.8943   -163.0393

Interpretación:

  • \(\hat{\beta}_1 = 120.8943\): Manteniendo las demás variables constantes, por cada semana adicional de gestación, ese espera un aumento en el peso al nacer de 120.8943 gramos.

  • \(\hat{\beta}_2 = -163.0393\): Manteniendo la edad gestacional constante, el hecho de ser mujer se asocia con una disminución esperada de 163.0393 gramos en el peso al nacer (en comparación con los hombres).

  1. Estime el peso esperado de una bebita con 40 semanas de gestación.
nuevo_dato <- data.frame(Edad = 40, Sexo = "Mujer")
peso_estimado <- predict(modelo_1.2, newdata = nuevo_dato, type = "response")
peso_estimado
##        1 
## 3062.451

Se estima que una bebé con 40 semanas de gestación tendría un peso aproximado de 3062.45 gramos.

Punto 2

Una muestra de personas de edad avanzada fue sometida a un examen psiquiátrico para determinar si estaban presentes los síntomas de senilidad. Al mismo tiempo, a cada persona se le aplicó la escala de inteligencia de Wechsler. A seguir se describe los puntajes obtenidos en la escala de Wechsler junto con la presencia (1) o ausencia (0) de síntomas de senilidad.

datos_senilidad <- data.frame(
  Escala = c(
    9, 13, 6, 8, 10, 4, 14, 8, 11, 7, 9, 7, 5, 14, 13, 16, 10, 12,
    11, 14, 15, 18, 7, 16, 9, 9, 11, 13, 15, 13, 10, 11, 6, 17, 14, 19,
    9, 11, 14, 10, 16, 16, 14, 13, 13, 9, 15, 10, 11, 12, 4, 14, 20, 10
  ),
  Senilidad = c(
    rep(1, 14), 
    rep(0, 40)
  )
)
  1. Ajuste modelos de respuesta binaria con funciones de enlace logit, probit y complemento log-log. Use el desvío y el criterio AIC para seleccionar el “mejor” modelo.
modelo_2.1 <- glm(
  Senilidad ~ Escala,
  family = binomial(link = "logit"),
  data = datos_senilidad
)

modelo_2.2 <- glm(
  Senilidad ~ Escala,
  family = binomial(link = "probit"),
  data = datos_senilidad
)

modelo_2.3 <- glm(
  Senilidad ~ Escala,
  family = binomial(link = "cloglog"),
  data = datos_senilidad
)

Desvio y criterios de informacion de los tres modelos

Modelo Deviance AIC
logit 51.01738 55.01738
probit 50.98358 54.98358
cloglog 51.31071 55.31071

El modelo con enlace probit se considera el “mejor” modelo, dado que presenta los criterios de ajuste más favorables

  1. Utilice las estadísticas de Wald y la razón de verosimilitudes para evaluar, al nivel de significancia aproximado del 5%, si la probabilidad de presentar síntomas de senilidad depende del puntaje en la escala de inteligencia de Wechsler.
summary(modelo_2.2)
## 
## Call:
## glm(formula = Senilidad ~ Escala, family = binomial(link = "probit"), 
##     data = datos_senilidad)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  1.38616    0.68528   2.023  0.04310 * 
## Escala      -0.18801    0.06301  -2.984  0.00284 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 61.806  on 53  degrees of freedom
## Residual deviance: 50.984  on 52  degrees of freedom
## AIC: 54.984
## 
## Number of Fisher Scoring iterations: 5

Interpretación

Una prueba de hipótesis del tipo \(H_0: \beta_1 = 0\) frente a \(H_1: \beta_1 \neq 0\) permite evaluar si la probabilidad de presentar síntomas de senilidad depende del puntaje obtenido en la prueba.

Como se observa en el resumen del modelo, el estadístico de Wald es \(-2.984\) y el valor-p correspondiente es \(0.00284\). Por lo tanto, existe evidencia estadísticamente significativa para rechazar la hipótesis nula de que \(\beta_1 = 0\). Esto indica que el puntaje en la escala de inteligencia de Wechsler tiene un efecto significativo sobre la probabilidad de presentar síntomas de senilidad.

La prueba de razón de verosimilitudes tiene como hipótesis nula \(H_0\) que el modelo nulo (\(\text{Modelo}_0\)) es igual al modelo completo (\(\text{Modelo}_1\)).

Dado que \[ \text{Modelo}_0: \beta_0 \quad \text{y} \quad \text{Modelo}_1: \beta_0 + \beta_1 X_{1i}, \] la hipótesis nula de esta prueba lleva a \[ H_0: \beta_1 = 0, \] Es decir, que la variable \(X_{1i}\) que es el puntaje en la escala de inteligencia de Wechsler no tiene efecto sobre la probabilidad de presentar senilidad.

library(lmtest)
modelo_2.2.0 <- glm(
  Senilidad ~ 1,
  family = binomial(link = "probit"),
  data = datos_senilidad
)

lrtest(modelo_2.2, modelo_2.2.0)
## Likelihood ratio test
## 
## Model 1: Senilidad ~ Escala
## Model 2: Senilidad ~ 1
##   #Df  LogLik Df  Chisq Pr(>Chisq)   
## 1   2 -25.492                        
## 2   1 -30.903 -1 10.823   0.001003 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretación:

Dado que \(p \approx 0.001 < 0.05\), rechazamos la hipótesis nula.

Por lo tanto, existe evidencia estadísticamente significativa de que el puntaje en la escala influye en la probabilidad de presentar senilidad

c)Realice el análisis de diagnóstico al modelo seleccionado. Comente.

library(performance)
check_model(modelo_2.2)

Se observa que los valores predichos por el modelo se encuentran bastante próximos a los valores realmente observados asi como para el grafico de cuantiles, lo cual sugiere un buen ajuste general. Sin embargo, en lo referente a los datos influyentes y tambien en las estimaciones de probabilidades pequeñas pareciera haber problemas.

Gráfico de envolvente (simulación):

El gráfico muestra un comportamiento satisfactorio, ya que los residuales se mantienen dentro de las bandas de confianza. Esto confirma que la distribución asumida por el modelo es correcta.

Distancia de Cook (valores influyentes):

plot(modelo_2.2, which = 4, main = "Distancia de Cook")

Al analizar la influencia de las observaciones, se destacan los datos número 7, 14 y 51, los cuales presentan una distancia de Cook que supera el criterio \(4/(n-p) \approx 0.07\). A continuación, se eliminará la observación con la distancia de Cook más alta para evaluar el impacto de este valor en la estimación del modelo.

ds1=datos_senilidad[-51, ]
modelo_2.2.1 <- glm(
  Senilidad ~  Escala,
  family = binomial(link = "probit"),
  data = ds1
)

summary(modelo_2.2.1)
## 
## Call:
## glm(formula = Senilidad ~ Escala, family = binomial(link = "probit"), 
##     data = ds1)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  1.84800    0.76658   2.411  0.01592 * 
## Escala      -0.22680    0.07071  -3.208  0.00134 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 61.200  on 52  degrees of freedom
## Residual deviance: 47.937  on 51  degrees of freedom
## AIC: 51.937
## 
## Number of Fisher Scoring iterations: 5

Dado que al eliminar la observación 51 la estimación de los parámetros cambia de manera notable, incluso alterando el signo de un coeficiente, se considerará el modelo sin esta observación. Posteriormente, se evaluará si la exclusión de las siguientes observaciones influyentes provoca cambios significativos en la estimación del modelo.

modelo_2.2.2 <- glm(
  Senilidad ~  Escala,
  family = binomial(link = "probit"),
  data = datos_senilidad[-c(7,14,51), ]
)

summary(modelo_2.2.2)
## 
## Call:
## glm(formula = Senilidad ~ Escala, family = binomial(link = "probit"), 
##     data = datos_senilidad[-c(7, 14, 51), ])
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.72129    0.97457   2.792 0.005234 ** 
## Escala      -0.33282    0.09792  -3.399 0.000677 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 55.651  on 50  degrees of freedom
## Residual deviance: 36.440  on 49  degrees of freedom
## AIC: 40.44
## 
## Number of Fisher Scoring iterations: 6

Dado que la exclusión de estas observaciones no genera cambios significativos en la estimación de los coeficientes, se considera que el modelo sin la observación 51 es el más adecuado.

  1. Interprete las estimaciones de los parámetros, excluyendo el intercepto.
modelo_2.2.1$coefficients
## (Intercept)      Escala 
##   1.8480049  -0.2268033

Interpretación: Dado que se está utilizando un modelo con enlace probit, la interpretación de las estimaciones no es completamente directa. No obstante, se puede afirmar que, a medida que aumenta el puntaje en la prueba, la probabilidad de presentar senilidad tiende a disminuir

  1. ¿Cuál es la probabilidad estimada de que una persona con 15 puntos en la escala de Wechsler presente síntomas de senilidad?
nuevo <- data.frame(Escala = 15)
prob_estimada <- predict(modelo_2.2.1, newdata = nuevo, type = "response")
prob_estimada
##          1 
## 0.06008692

Se estima que la probabilidad de que una persona con 15 puntos en la escala de Wechsler presente síntomas de senilidad es 0.06008692

Punto 3

Estos datos son producidos por un experimento donde lotes de alrededor de cincuenta insectos fueron rociados con diferentes concentraciones del plaguicida conocido como rotenona. Los insectos aparentemente muertos, moribundos o tan gravemente afectados como para ser incapaces de caminar más de unos pocos pasos se clasificaron como muertos.

datos_insectos <- data.frame(
  Dosis = c(0.41, 0.58, 0.71, 0.89, 1.01),
  Expuestos = c(50, 48, 46, 49, 50),
  Muertos = c(6, 16, 24, 42, 44)
)
  1. Realice un gráfico de la tasa de mortalidad de los insectos como función de la concentración del plaguicida. Comente.

La tasa de mortalidad de los insectos aumenta conforme crece la concentración del plaguicida, mostrando una relación directa y clara entre dosis y efecto letal. En concentraciones altas, la mortalidad supera el 85%, indicando un fuerte impacto del plaguicida

  1. Ajuste modelos de respuesta binaria con funciones de enlace logit, probit y complemento log-log. Compare gráficamente los tres modelos ajustados. Comente.
respuesta <- cbind(datos_insectos$Muertos, datos_insectos$Expuestos - datos_insectos$Muertos)

modelo_logit <- glm(respuesta ~ Dosis, family = binomial(link = "logit"), data = datos_insectos)
modelo_probit <- glm(respuesta ~ Dosis, family = binomial(link = "probit"), data = datos_insectos)
modelo_cloglog <- glm(respuesta ~ Dosis, family = binomial(link = "cloglog"), data = datos_insectos)

Las estimaciones usando los tres enlaces son bastante similares, no hay diferencias muy notorias.

  1. Compare los tres modelos ajustados usando los valores del desvío y del AIC. Seleccione el ``mejor’’ modelo.

Desvio y criterios de informacion de los tres modelos

Modelo Deviance AIC
logit 1.348060 24.56846
probit 1.638635 24.85903
cloglog 4.035474 27.25587

El modelo con enlace logit se considera el “mejor” modelo, dado que presenta los criterios de ajuste más favorables

  1. Use las estadísticas de Wald y razón de verosimilitudes para evaluar, al nivel de significancia aproximado de 5 %, si la concentración del plaguicida afecta la probabilidad que los insectos mueran.

Interpretación:

summary(modelo_logit)
## 
## Call:
## glm(formula = respuesta ~ Dosis, family = binomial(link = "logit"), 
##     data = datos_insectos)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -4.8387     0.6377  -7.588 3.24e-14 ***
## Dosis         7.0681     0.8826   8.009 1.16e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 96.6881  on 4  degrees of freedom
## Residual deviance:  1.3481  on 3  degrees of freedom
## AIC: 24.568
## 
## Number of Fisher Scoring iterations: 4

Una prueba de hipótesis del tipo \[ H_0: \beta_1 = 0 \qquad \text{vs.} \qquad H_1: \beta_1 \neq 0 \] Permite evaluar si la concentración del plaguicida influye en la probabilidad de que los insectos mueran.

En el resumen del modelo, el estadístico de Wald es \(8.009\) y el valor \(p\) correspondiente es \(1.16 \times 10^{-15}\). Por lo tanto, existe evidencia estadísticamente significativa para rechazar la hipótesis nula \((\beta_1 = 0)\). Esto indica que la concentración del plaguicida sí afecta la probabilidad de mortalidad de los insectos.

La prueba de razón de verosimilitudes contrasta como hipótesis nula que el modelo nulo (\(\text{Modelo}_0\)) es equivalente al modelo completo (\(\text{Modelo}_1\)). Dado que \[ \text{Modelo}_0: \beta_0 \qquad \text{Modelo}_1: \beta_0 + \beta_1 X_i, \] La hipótesis nula se puede expresarse como \[ H_0: \beta_1 = 0, \] Es decir, que la concentración del plaguicida no tiene efecto sobre la probabilidad de muerte de los insectos.

library(lmtest)
modelo_logit.0 <- glm(
  respuesta ~ 1,
  family = binomial(link = "logit"),
  data = datos_insectos
)

lrtest(modelo_logit, modelo_logit.0)
## Likelihood ratio test
## 
## Model 1: respuesta ~ Dosis
## Model 2: respuesta ~ 1
##   #Df  LogLik Df Chisq Pr(>Chisq)    
## 1   2 -10.284                        
## 2   1 -57.954 -1 95.34  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretación: Dado que \(p \ll 0.05\), se rechaza la hipótesis nula. En consecuencia, existe evidencia estadísticamente significativa de que la concentración del plaguicida afecta la probabilidad de que los insectos mueran.

  1. Realice el análisis de diagnóstico al modelo. Comente.
library(performance)
check_model(modelo_logit)

Se observa que el gráfico de homogeneidad no presenta el comportamiento esperado, ya que los residuos no siguen una línea recta. En cuanto a los valores influyentes, se aprecian problemas. El gráfico de cuantiles no muestra preocupaciones relevantes; sin embargo, debido a la pequeña cantidad de datos disponibles, este tipo de validaciones tiene una capacidad limitada.

Gráfico de envolvente (simulación):

El gráfico muestra un comportamiento satisfactorio, ya que los residuales se mantienen dentro de las bandas de confianza.

Distancia de Cook (valores influyentes):

plot(modelo_logit, which = 4, main = "Distancia de Cook")

Al analizar la influencia de las observaciones, es importante considerar que este tipo de verificaciones puede resultar limitado o poco fiable debido a la reducida cantidad de datos disponibles.

  1. Estime la dosis letal 50 % y calcule para ella un intervalo de confianza de aproximadamente 95 %.

La relacion entre la media y el predictor lineal se expresa como \[ \operatorname{logit}(\mu) = \beta_0 + \beta_1 X_i, \] al evaluar en \(\mu = 0.5\), dado que \(\operatorname{logit}(0.5) = 0\)se obtiene \[ 0 = \beta_0 + \beta_1 X_i. \] Por lo tanto, la dosis letal al 50% (DL50) viene dada por \[ X_i = -\frac{\beta_0}{\beta_1}. \]

b0 <- coef(modelo_logit)[1]
b1 <- coef(modelo_logit)[2]

DL50 <- -b0 / b1
DL50
## (Intercept) 
##   0.6845816
# Usando la libreria MASS
library(MASS)
IC <- dose.p(modelo_logit, p=0.5)  
IC
##               Dose         SE
## p = 0.5: 0.6845816 0.02271034

La estimación puntual de la DL50 fue

\(\hat{\theta} = 0.6846\) con un error estándar de \(SE = 0.0227\). El intervalo de confianza del 95% es: \[ 0.6846 \pm 1.96 \times 0.0227 = (0.6401,\; 0.7291). \]

Punto 4

Se desarrolla un experimento donde ratones hembras fueron alimentados con dosis extremadamente bajas de un conocido carcinógeno denominado 2-acetilaminofluoreno. A seguir se describe el número de ratones expuestos a cada dosis del carcinógeno y el número de ellos que desarrollaron cáncer.

datos_higado <- data.frame(
  Dosis = c(0.00, 0.30, 0.35, 0.45, 0.60, 0.75, 1.00, 1.50),
  Expuestos = c(555, 2014, 1102, 550, 441, 382, 213, 211),
  Enfermos = c(6, 34, 20, 15, 13, 17, 19, 24)
)

datos_vejiga <- data.frame(
  Dosis = c(0.00, 0.30, 0.35, 0.45, 0.60, 0.75, 1.00, 1.50),
  Expuestos = c(101, 443, 200, 103, 66, 75, 31, 11),
  Enfermos = c(1, 5, 0, 2, 2, 12, 21, 11)
)
  1. Analice los datos para cada tipo de cáncer separadamente. Ajuste modelos de respuesta binaria (enlace logit) y respuesta Poisson (enlace log) según sea el caso. Justifique la elección del modelo.
resp_higado4 <- cbind(datos_higado$Enfermos, datos_higado$Expuestos - datos_higado$Enfermos)

resp_vejiga4 <- cbind(datos_vejiga$Enfermos, datos_vejiga$Expuestos - datos_vejiga$Enfermos)

Modelo para cancer de higado

modelo_higado_bin <- glm(
  resp_higado4 ~ Dosis,
  family = binomial(link = "logit"),
  data = datos_higado
)
modelo_higado_pois <- glm(
  Enfermos ~ Dosis + offset(log(Expuestos)),
  family = poisson(link = "log"),
  data = datos_higado
)

Modelo para cancer de vejiga

modelo_vejiga_bin <- glm(
  resp_vejiga4 ~ Dosis,
  family = binomial(link = "logit"),
  data = datos_vejiga
)
modelo_vejiga_poi <- glm(
  Enfermos ~ Dosis + offset(log(Expuestos)),
  family = poisson(link = "log"),
  data = datos_vejiga
)

El análisis se realiza de manera separada para cada tipo de cáncer, ya que los mecanismos biológicos y la relación con la dosis pueden diferir entre el cáncer de hígado y el cáncer de vejiga.

La elección del modelo ya sea de respuesta binomial o poisson depende del objetivo del análisis. Si el interés se centra en estimar la probabilidad de desarrollar cáncer en función de la dosis, se utiliza un modelo binomial. En cambio, si el objetivo es modelar la tasa de ocurrencia del cáncer, es decir, el número de casos por unidad de exposición, se emplea un modelo con respuesta poisson.

  1. Use las estadísticas de Wald y razón de verosimilitudes para evaluar, al nivel de significancia aproximado de 5%, si la dosis del carcinógeno afecta la probabilidad/incidencia de cáncer en los ratones.

Modelo para cancer de higado

summary(modelo_higado_bin)
## 
## Call:
## glm(formula = resp_higado4 ~ Dosis, family = binomial(link = "logit"), 
##     data = datos_higado)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -4.5083     0.1478 -30.501   <2e-16 ***
## Dosis         1.7624     0.1864   9.457   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 77.6029  on 7  degrees of freedom
## Residual deviance:  4.5217  on 6  degrees of freedom
## AIC: 45.494
## 
## Number of Fisher Scoring iterations: 4
summary(modelo_higado_pois)
## 
## Call:
## glm(formula = Enfermos ~ Dosis + offset(log(Expuestos)), family = poisson(link = "log"), 
##     data = datos_higado)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -4.4919     0.1445 -31.096   <2e-16 ***
## Dosis         1.6634     0.1767   9.414   <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: 74.4827  on 7  degrees of freedom
## Residual deviance:  4.6882  on 6  degrees of freedom
## AIC: 46.024
## 
## Number of Fisher Scoring iterations: 4

Tanto en el modelo binomial (probabilidad) como en el modelo de Poisson (incidencia), el coeficiente asociado a la dosis del carcinógeno resulta , el valor \(z\) para la dosis es aproximadamente 9.4 en ambos modelos, con \[ p < 2 \times 10^{-16} \ll 0.05, \] Lo que indica un efecto significativo de la dosis.

library(lmtest)
modelo_higado_bin.0 <- glm(
  resp_higado4 ~ 1,
  family = binomial(link = "logit"),
  data = datos_higado
)

modelo_higado_pois.0 <- glm(
  Enfermos ~  offset(log(Expuestos)),
  family = poisson(link = "log"),
  data = datos_higado
)

lrtest(modelo_higado_bin, modelo_higado_bin.0)
## Likelihood ratio test
## 
## Model 1: resp_higado4 ~ Dosis
## Model 2: resp_higado4 ~ 1
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   2 -20.747                         
## 2   1 -57.288 -1 73.081  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(modelo_higado_pois, modelo_higado_pois.0)
## Likelihood ratio test
## 
## Model 1: Enfermos ~ Dosis + offset(log(Expuestos))
## Model 2: Enfermos ~ offset(log(Expuestos))
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   2 -21.012                         
## 2   1 -55.909 -1 69.794  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La prueba muestra que la inclusión de la dosis del carcinógeno es significativa tanto del modelo binomial (probabilidad de cáncer) como del modelo Poisson con offset (incidencia de cáncer).

En ambos casos, el estadístico de prueba es muy alto (\(\chi^2 \approx 73\) para el modelo binomial y \(\chi^2 \approx 69.8\) para el modelo Poisson), con , es decir, mucho menor al nivel de significancia de 0.05.

Por lo tanto, existe evidencia estadísticamente significativa de que la dosis del carcinógeno afecta tanto la probabilidad como la incidencia de cáncer de hígado en los ratones. Esto confirma que la variable dosis es un factor determinante en la aparición del cáncer.

Modelo para cancer de vejiga

summary(modelo_vejiga_bin)
## 
## Call:
## glm(formula = resp_vejiga4 ~ Dosis, family = binomial(link = "logit"), 
##     data = datos_vejiga)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -7.4324     0.5831  -12.75   <2e-16 ***
## Dosis         7.8752     0.7821   10.07   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 214.94  on 7  degrees of freedom
## Residual deviance:  11.45  on 6  degrees of freedom
## AIC: 34.019
## 
## Number of Fisher Scoring iterations: 5
summary(modelo_vejiga_poi)
## 
## Call:
## glm(formula = Enfermos ~ Dosis + offset(log(Expuestos)), family = poisson(link = "log"), 
##     data = datos_vejiga)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.2178     0.2942  -17.74   <2e-16 ***
## Dosis         3.8602     0.2811   13.73   <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: 174.256  on 7  degrees of freedom
## Residual deviance:  32.564  on 6  degrees of freedom
## AIC: 60.75
## 
## Number of Fisher Scoring iterations: 5

Tanto en el modelo binomial (probabilidad) como en el modelo de Poisson (incidencia), el coeficiente asociado a la dosis del carcinógeno resulta , el valor \(z\) para la dosis es aproximadamente 10.07 y 13.73 respectivamente. \[ p < 2 \times 10^{-16} \ll 0.05, \] Lo que indica un efecto significativo de la dosis.

library(lmtest)
modelo_vejiga_bin.0 <- glm(
  resp_vejiga4 ~ 1,
  family = binomial(link = "logit"),
  data = datos_vejiga
)

modelo_vejiga_poi.0 <- glm(
  Enfermos ~  offset(log(Expuestos)),
  family = poisson(link = "log"),
  data = datos_vejiga
)

lrtest(modelo_vejiga_bin, modelo_vejiga_bin.0)
## Likelihood ratio test
## 
## Model 1: resp_vejiga4 ~ Dosis
## Model 2: resp_vejiga4 ~ 1
##   #Df  LogLik Df Chisq Pr(>Chisq)    
## 1   2  -15.01                        
## 2   1 -116.76 -1 203.5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lrtest(modelo_vejiga_poi, modelo_vejiga_poi.0)
## Likelihood ratio test
## 
## Model 1: Enfermos ~ Dosis + offset(log(Expuestos))
## Model 2: Enfermos ~ offset(log(Expuestos))
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   2 -28.375                         
## 2   1 -99.221 -1 141.69  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretacion La prueba muestra que la inclusión de la dosis del carcinógeno en los modelos mejora significativamente el ajuste tanto del modelo binomial (probabilidad de cáncer) como del modelo Poisson con offset (incidencia de cáncer).

En ambos casos, el estadístico de prueba es muy alto (\(\chi^2 \approx 203.5\) para el modelo binomial y \(\chi^2 \approx 141.69\) para el modelo Poisson), con , es decir, mucho menor al nivel de significancia de 0.05.

Por lo tanto, existe evidencia estadísticamente significativa de que la dosis del carcinógeno afecta tanto la probabilidad como la incidencia de cáncer de vejiga en los ratones. Esto confirma que la variable dosis es un factor determinante en la aparición del cáncer.

  1. Realice el análisis de diagnóstico a los modelos. Comente.
library(performance)
check_model(modelo_higado_bin)

Gráfico de envolvente (simulación):

library(performance)
check_model(modelo_higado_pois)

Gráfico de envolvente (simulación):

library(performance)
check_model(modelo_vejiga_bin)

Gráfico de envolvente (simulación):

library(performance)
check_model(modelo_vejiga_poi)

Gráfico de envolvente (simulación):

Interpretación:

El diagnóstico de los modelos pone en evidencia la presencia de un dato influyente, correspondiente a la observación número 8. Asimismo, en el gráfico de homocedasticidad se observan algunas particularidades. En particular, para los modelos ajustados con respuesta Poisson en el caso del cáncer de vejiga, se aprecia un posible problema en la dispersión, ya que los valores se alejan considerablemente de las líneas de referencia de los valores predichos y observados. Este comportamiento también se ve reflejado en la relación entre la devianza residual y los grados de libertad para este caso específico. No obstante, los gráficos envelope no muestran problemas graves, lo cual puede deberse, en parte, al reducido número de observaciones disponibles.

  1. Interprete las estimaciones de los parámetros excepto los interceptos.

Un aumento de 1 unidades en las dosis, aumenta los chances promedio de tener cancer de higado y vejiga en 19% y 119% respectivamente.

Punto 5

Los siguientes datos corresponden a las fallas de piezas de equipos electrónicos que funcionan en dos modos. Para cada pieza electrónica se registra el tiempo de funcionamiento en el modo 1 y el modo 2 así como el número total de fallas.

datos_fallas <- data.frame(
  Tiempo_Modo1 = c(33.3, 52.2, 64.7, 137.0, 125.9, 116.3, 131.7, 85.0, 91.9),
  Tiempo_Modo2 = c(25.3, 14.4, 32.5, 20.5, 97.6, 53.6, 56.6, 87.3, 47.8),
  Numero_Fallas = c(15, 9, 14, 24, 27, 27, 23, 18, 22)
)
  1. Ajuste modelos de respuesta Poisson con función de enlace logarítmica, identidad y raíz cuadrada, donde el número de fallas es la variable respuesta y los tiempos de funcionamiento en el modo 1 y 2 son las variables de explicación. Use el desvío y el criterio AIC para seleccionar el ``mejor’’ modelo.
modelo_fallas_log <- glm(
  Numero_Fallas ~  Tiempo_Modo1 + Tiempo_Modo2,
  family = poisson(link = "log"),
  data = datos_fallas
)
modelo_fallas_ide <- glm(
  Numero_Fallas ~  Tiempo_Modo1 + Tiempo_Modo2,
  family = poisson(link = "identity"),
  data = datos_fallas
)
modelo_fallas_raiz <- glm(
  Numero_Fallas ~  Tiempo_Modo1 + Tiempo_Modo2,
  family = poisson(link = "sqrt"),
  data = datos_fallas
)

Desvio y criterios de informacion de los tres modelos

Modelo Deviance AIC
log 4.003304 53.06006
identidad 4.197149 53.25390
raiz 4.034107 53.09086

El modelo con enlace log se considera el “mejor” modelo, dado que presenta los criterios de ajuste más favorables

  1. Use las estadísticas de Wald y razón de verosimilitudes para evaluar, al nivel de significancia aproximado de \(5\%\), si el número de fallas de las piezas electrónicas depende del tiempo de funcionamiento en el modo 1.
summary(modelo_fallas_log)
## 
## Call:
## glm(formula = Numero_Fallas ~ Tiempo_Modo1 + Tiempo_Modo2, family = poisson(link = "log"), 
##     data = datos_fallas)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.175168   0.255456   8.515  < 2e-16 ***
## Tiempo_Modo1 0.007015   0.002429   2.888  0.00387 ** 
## Tiempo_Modo2 0.002549   0.002835   0.899  0.36852    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 16.9964  on 8  degrees of freedom
## Residual deviance:  4.0033  on 6  degrees of freedom
## AIC: 53.06
## 
## Number of Fisher Scoring iterations: 4
modelo_fallas_log.0 <- glm(
  Numero_Fallas ~  Tiempo_Modo2,
  family = poisson(link = "log"),
  data = datos_fallas
)

lrtest(modelo_fallas_log,modelo_fallas_log.0)
## Likelihood ratio test
## 
## Model 1: Numero_Fallas ~ Tiempo_Modo1 + Tiempo_Modo2
## Model 2: Numero_Fallas ~ Tiempo_Modo2
##   #Df  LogLik Df  Chisq Pr(>Chisq)   
## 1   3 -23.530                        
## 2   2 -27.706 -1 8.3528   0.003851 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A partir del resumen del modelo, el valor p obtenido es 0.00387, el cual es menor que el nivel de significancia de 0.05, se concluye que existe evidencia estadísticamente significativa para afirmar que el tiempo en modo 1 tiene un efecto sobre el número de fallas de las piezas electrónicas. Asimismo, mediante la prueba de razón de verosimilitudes se obtiene un valor p de 0.003, igualmente inferior al nivel de significancia considerado, lo que refuerza esta conclusión. Por lo tanto, por ambos métodos se confirma que el coeficiente asociado al tiempo en modo 1 es estadísticamente significativo en el modelo, lo que indica que esta variable explica de manera significativa la variabilidad en el número de fallas de las piezas electrónicas.

  1. Use las estadísticas de Wald y razón de verosimilitudes para evaluar, al nivel de significancia aproximado de \(5\%\), si el número de fallas de las piezas electrónicas depende del tiempo de funcionamiento en el modo 2.
modelo_fallas_log.00 <- glm(
  Numero_Fallas ~  Tiempo_Modo1,
  family = poisson(link = "log"),
  data = datos_fallas
)

lrtest(modelo_fallas_log,modelo_fallas_log.00)
## Likelihood ratio test
## 
## Model 1: Numero_Fallas ~ Tiempo_Modo1 + Tiempo_Modo2
## Model 2: Numero_Fallas ~ Tiempo_Modo1
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   3 -23.530                     
## 2   2 -23.932 -1 0.8045     0.3697

Para analizar el coeficiente asociado al tiempo en modo 2, a partir del resumen del modelo se obtiene que este no es estadísticamente significativo, ya que su valor p es mayor que el nivel de significancia. De igual manera, al realizar la prueba de razón de verosimilitudes se obtiene un valor p de 0.36, el cual también es superior al nivel de significancia establecido, por lo que no se rechaza la hipótesis nula. En consecuencia, no se encuentra evidencia estadística suficiente para afirmar que \(\beta_2\) sea distinto de cero, y por lo tanto, el tiempo en modo 2 no explica de manera significativa el número de fallas de los aparatos electrónicos.

  1. Realice el análisis de diagnóstico al modelo. Comente.
library(performance)
check_model(modelo_fallas_log)

Al realizar el diagnóstico del modelo se observa que las predicciones caen dentro de los intervalos de los valores observados. No obstante, hay problemas con la dispersión: las curvas de distribución observada y simulada no siguen la misma forma, lo que indica discrepancias en la varianza entre los datos y lo que el modelo supone. En el gráfico de homogeneidad de varianza la tendencia no es completamente horizontal (hay desviaciones), aunque dichas desviaciones no son extremadamente grandes. Además, se identifican varios valores que podrían ser influyentes.

Gráfico de envolvente (simulación):

Distancia de Cook (valores influyentes):

plot(modelo_fallas_log, which = 4, main = "Distancia de Cook")

  1. Interprete las estimaciones de los parámetros excepto el intercepto.

Por cada aumento de una unidad en el tiempo de funcionamiento del modo 1, el numero de fallas promedio aumenta en 0.77%.

Punto 6

Un investigador contó el número de lesiones causadas por cierto virus después de su exposición a los rayos X durante varios minutos.

datos_virus <- data.frame(
  Tiempo = c(0, 15, 30, 45, 60),
  Lesiones = c(271, 108, 59, 29, 12)
)
  1. Realice un gráfico del número de lesiones como función del tiempo de exposición. Comente.
plot(datos_virus$Tiempo, datos_virus$Lesiones,
     type = "b", pch = 19,
     xlab = "Tiempo de exposición",
     ylab = "Número de lesiones",
     main = "Lesiones en función del tiempo de exposición")

  1. Ajuste modelos de respuesta Poisson con función de enlace logarítmica, identidad y raíz cuadrada. Use el desvío y el criterio AIC para seleccionar el ``mejor’’ modelo.
modelo_virus_log <- glm(
  Lesiones ~  Tiempo,
  family = poisson(link = "log"),
  data = datos_virus
)
modelo_virus_ide <- glm(
  Lesiones ~  Tiempo,
  family = poisson(link = "identity"),
  data = datos_virus
)
modelo_virus_raiz <- glm(
  Lesiones ~  Tiempo,
  family = poisson(link = "sqrt"),
  data = datos_virus
)

Desvio y criterios de informacion de los tres modelos

Modelo Deviance AIC
log 2.290581 35.71857
identidad 74.184372 107.61236
raiz 25.946845 59.37484

El modelo con enlace log se considera el “mejor” modelo, dado que presenta los criterios de ajuste más favorables

  1. Use las estadísticas de Wald y razón de verosimilitudes para evaluar, al nivel de significancia aproximado de \(5\%\), si el número de lesiones depende del tiempo de exposición.
summary(modelo_virus_log)
## 
## Call:
## glm(formula = Lesiones ~ Tiempo, family = poisson(link = "log"), 
##     data = datos_virus)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.571340   0.056714   98.24   <2e-16 ***
## Tiempo      -0.051326   0.002972  -17.27   <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: 413.1314  on 4  degrees of freedom
## Residual deviance:   2.2906  on 3  degrees of freedom
## AIC: 35.719
## 
## Number of Fisher Scoring iterations: 3
modelo_virus_log.0 <- glm(
  Lesiones ~  1,
  family = poisson(link = "log"),
  data = datos_virus
)

lrtest(modelo_virus_log,modelo_virus_log.0)
## Likelihood ratio test
## 
## Model 1: Lesiones ~ Tiempo
## Model 2: Lesiones ~ 1
##   #Df   LogLik Df  Chisq Pr(>Chisq)    
## 1   2  -15.859                         
## 2   1 -221.280 -1 410.84  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A partir del resumen del modelo, se observa que la variable tiempo tiene un p-valor inferior al nivel de significancia \(\alpha=0.05\), por lo que existe evidencia estadísticamente significativa para rechazar la hipótesis nula de que \(\beta_1 =0\). Esto indica que el tiempo es una variable explicativa relevante del número de lesiones. Adicionalmente, la prueba de razón de verosimilitudes confirma este resultado, reforzando la evidencia de un efecto significativo del tiempo sobre la respuesta.

  1. Realice el análisis de diagnóstico al modelo. Comente.
library(performance)
check_model(modelo_virus_log)

Gráfico de envolvente (simulación):

Distancia de Cook (valores influyentes):

plot(modelo_virus_log, which = 4, main = "Distancia de Cook")

Al realizar el diagnóstico del modelo, se observa que la mayoría de los valores observados se encuentran dentro de los intervalos predichos por el modelo, lo cual es un resultado favorable. No obstante, se evidencian problemas en la gráfica de dispersión, ya que las dos curvas no siguen un patrón similar. Asimismo, en la presencia de observaciones influyentes, por otro lado, el gráfico de envolvente muestra un comportamiento adecuado.

Finalmente, es importante considerar que el tamaño muestral es muy reducido, lo cual puede afectar la interpretación de los resultados y limitar la validez del modelo, debido a la escasa cantidad de datos disponibles.

e)Interprete las estimaciones de los parámetros excepto el intercepto.

  1. Estime el número esperado de lesiones causadas por el virus después de su exposición a los 10 días.
nuevo <- data.frame(Tiempo = 10)
prob_estimada <- predict(modelo_virus_log, newdata = nuevo, type = "response")
prob_estimada
##       1 
## 157.289

Punto 7

Estos datos muestran la incidencia de cáncer de piel por rangos de edad en las mujeres de las ciudades de St Paul (Minneapolis-Minnesota) y Fort Worth (Dallas-Texas). Se espera que la incidencia de este tipo de cáncer sea mayor en Texas que en Minnesota ya que la exposición al sol es mayor en Texas.

datos_cancer <- data.frame(
  Ciudad = factor(c(
    rep("St Paul", 7),
    rep("Forth Worth", 7)
  )),

  Edad = factor(rep(c(
    "15-24", "25-34", "35-44", "45-54", "55-64", "65-74", "75+"
  ), 2)),

  Casos = c(
    1, 16, 30, 71, 102, 130, 40,
    4, 38, 119, 221, 259, 310, 65
  ),

  Poblacion = c(
    172675, 123065, 96216, 92051, 72159, 54722, 8328,
    181343, 146207, 121374, 111353, 83004, 55932, 7583
  )
)
  1. Ajuste un modelo de respuesta Poisson con función de enlace logarítmica y las covariables cualitativas ciudad y edad. Comente.
modelo_cancer_log <- glm(
  Casos ~ Edad + Ciudad + offset(log(Poblacion)),
  family = poisson(link = "log"),
  data = datos_cancer
)

summary(modelo_cancer_log)
## 
## Call:
## glm(formula = Casos ~ Edad + Ciudad + offset(log(Poblacion)), 
##     family = poisson(link = "log"), data = datos_cancer)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -10.83938    0.44754 -24.220  < 2e-16 ***
## Edad25-34       2.62899    0.46746   5.624 1.87e-08 ***
## Edad35-44       3.84558    0.45466   8.458  < 2e-16 ***
## Edad45-54       4.59381    0.45103  10.185  < 2e-16 ***
## Edad55-64       5.08638    0.45030  11.296  < 2e-16 ***
## Edad65-74       5.64569    0.44975  12.553  < 2e-16 ***
## Edad75+         6.17568    0.45774  13.492  < 2e-16 ***
## CiudadSt Paul  -0.85269    0.05962 -14.302  < 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: 2169.8429  on 13  degrees of freedom
## Residual deviance:    5.2089  on  6  degrees of freedom
## AIC: 101.47
## 
## Number of Fisher Scoring iterations: 4
  1. Use las estadísticas de Wald y razón de verosimilitudes para evaluar, al nivel de significancia aproximado de \(5\%\), si el cáncer de piel depende de la ciudad.
modelo_cancer_log <- glm(
  Casos ~ Edad + Ciudad + offset(log(Poblacion)),
  family = poisson(link = "log"),
  data = datos_cancer
)

modelo_cancer_log.0 <- glm(
  Casos ~ Edad + offset(log(Poblacion)),
  family = poisson(link = "log"),
  data = datos_cancer
)

lrtest(modelo_cancer_log,modelo_cancer_log.0)
## Likelihood ratio test
## 
## Model 1: Casos ~ Edad + Ciudad + offset(log(Poblacion))
## Model 2: Casos ~ Edad + offset(log(Poblacion))
##   #Df   LogLik Df  Chisq Pr(>Chisq)    
## 1   8  -42.733                         
## 2   7 -155.992 -1 226.52  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El test de Wald como la prueba de razon de verosimilitudes confirman que el parametro de la ciudad es significativo, es decir cancer de piel depende de la ciudad

  1. Realice el análisis de diagnóstico. Comente.
library(performance)
check_model(modelo_cancer_log,)

Gráfico de envolvente (simulación):

Distancia de Cook (valores influyentes):

plot(modelo_cancer_log,, which = 4, main = "Distancia de Cook")

Las graficas no muestran problemas muy significativos del modelo.

  1. ¿La incidencia de cáncer de piel es mayor en Fort Worth que en St Paul? Justifique.

Ya que el coeficiente de St.Paul es negativo, indica que a comparación de Forth Worth la tasas de cáncer de piel en promedio es menor ahí.

  1. Estime el número esperado de casos de cáncer de piel en personas de más de \(75\) años de edad en la ciudad de St Paul si la población en este rango de edad crece hasta 10000 personas.
nuevo=data.frame(Ciudad="St Paul",Edad="75+",Poblacion=10000)
predict(modelo_cancer_log,newdata = nuevo,type="response")
##        1 
## 40.20308

Punto 8

Estos datos describen una etapa en el proceso de fabricación de semiconductores. Se cree que cuatro factores influyen en la resistencia de la placa de los semiconductores. Por lo tanto, se realizó un experimento factorial completo con dos niveles de cada factor. La experiencia previa sugiere que la resistencia tiene una distribución estrictamente positiva y sesgada a la derecha. Los datos están disponibles en el objeto wafer del paquete faraway de R.

  1. Ajuste modelos de respuesta Gamma con funciones de enlace identidad, logarítmica y recíproca, y los cuatro factores como covariables. Use el criterio AIC para seleccionar el ``mejor’’ modelo.
library(faraway)
data(wafer)
head(wafer)
##   x1 x2 x3 x4 resist
## 1  -  -  -  -  193.4
## 2  +  -  -  -  247.6
## 3  -  +  -  -  168.2
## 4  +  +  -  -  205.0
## 5  -  -  +  -  303.4
## 6  +  -  +  -  339.9
# Modelo con enlace identidad
mod_id <- glm(resist ~ x1 + x2 + x3 + x4,
              family = Gamma(link = "identity"),
              data = wafer)

# Modelo con enlace logarítmico
mod_log <- glm(resist ~ x1 + x2 + x3 + x4,
               family = Gamma(link = "log"),
               data = wafer)

# Modelo con enlace recíproca
mod_inv <- glm(resist ~ x1 + x2 + x3 + x4,
               family = Gamma(link = "inverse"),
               data = wafer)

AIC(mod_id, mod_log, mod_inv)
##         df      AIC
## mod_id   6 154.8398
## mod_log  6 152.9078
## mod_inv  6 151.7337

El modelo con enlace reciproca se considera el “mejor” modelo, dado que presenta el criterio de ajuste más favorables

  1. Evalúe, al nivel de significancia aproximado de \(5\%\), cuál(es) factor(es) afecta(n) la resistencia de las placas de los semiconductores.
summary(mod_inv)
## 
## Call:
## glm(formula = resist ~ x1 + x2 + x3 + x4, family = Gamma(link = "inverse"), 
##     data = wafer)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0043603  0.0002471  17.644 2.04e-09 ***
## x1+         -0.0004756  0.0002180  -2.181 0.051752 .  
## x2+          0.0013445  0.0002272   5.918 0.000101 ***
## x3+         -0.0008148  0.0002208  -3.691 0.003560 ** 
## x4+          0.0002662  0.0002170   1.227 0.245574    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.01022894)
## 
##     Null deviance: 0.69784  on 15  degrees of freedom
## Residual deviance: 0.11540  on 11  degrees of freedom
## AIC: 151.73
## 
## Number of Fisher Scoring iterations: 4

Las variables significativas son x2 y x3 ya que el valor p es menor al nivel de significancia

  1. Realice el análisis de diagnóstico. Comente.
library(performance)
check_model(mod_inv)

Gráfico de envolvente (simulación):

Distancia de Cook (valores influyentes):

plot(mod_inv, which = 4, main = "Distancia de Cook")

No se evidencias problemas draticos en la especificacion del modelos