Taller Global Modelos Lineales Generalizados
Sebastian Romero
2025-11-27
Punto 1
Se cargan los datos y se verifica que el tipo de variables esten bien.
library(performance)
library(glmtoolbox)
library(car)## Cargando paquete requerido: carDatalibrary(MASS)
library(faraway)##
## Adjuntando el paquete: 'faraway'## The following objects are masked from 'package:car':
##
## logit, vif##Punto 1- ----------------------
datos1 <- read.table(header = TRUE, text = "
Sexo Edad Peso
Mujer 40 3317
Mujer 36 2729
Mujer 40 2935
Mujer 38 2754
Mujer 42 3210
Mujer 39 2817
Mujer 40 3126
Mujer 37 2539
Mujer 36 2412
Mujer 38 2991
Mujer 39 2875
Mujer 40 3231
Hombre 40 2968
Hombre 38 2795
Hombre 40 3163
Hombre 35 2925
Hombre 36 2625
Hombre 37 2847
Hombre 41 3292
Hombre 40 3473
Hombre 37 2628
Hombre 38 3176
Hombre 40 3421
Hombre 38 2975
")
head(datos1)## Sexo Edad Peso
## 1 Mujer 40 3317
## 2 Mujer 36 2729
## 3 Mujer 40 2935
## 4 Mujer 38 2754
## 5 Mujer 42 3210
## 6 Mujer 39 2817str(datos1)## 'data.frame': 24 obs. of 3 variables:
## $ Sexo: chr "Mujer" "Mujer" "Mujer" "Mujer" ...
## $ Edad: int 40 36 40 38 42 39 40 37 36 38 ...
## $ Peso: int 3317 2729 2935 2754 3210 2817 3126 2539 2412 2991 ...#Cambiar el tipo de la variable de Sexo
datos1$Sexo <- as.factor(datos1$Sexo)a)
Ajuste un modelo de respuesta normal con funcion de enlace identidad y las covariables edad gestacional, género y la interacci´on entre estas.
fit1=lm(Peso ~Sexo*Edad,data = datos1)
summary(fit1)##
## Call:
## lm(formula = Peso ~ Sexo * Edad, data = datos1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -246.69 -138.11 -39.13 176.57 274.28
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1268.67 1114.64 -1.138 0.268492
## SexoMujer -872.99 1611.33 -0.542 0.593952
## Edad 111.98 29.05 3.855 0.000986 ***
## SexoMujer:Edad 18.42 41.76 0.441 0.663893
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 180.6 on 20 degrees of freedom
## Multiple R-squared: 0.6435, Adjusted R-squared: 0.59
## F-statistic: 12.03 on 3 and 20 DF, p-value: 0.000101b)
Evalúe, al nivel de significancia de 5%, si el efecto de la edad gestacional sobre el peso al nacer es el mismo para niños y niñas.
Para evaluar esta hipótesis basta ver el p-valor respecto a la interacción, se aprecia que con una significania del 5%, la interacción no es significativa.
c)
Realice el análisis de diagnóstico. Comente.
check_model(fit1)
set.seed(123)
envelope(fit1, rep=1000)
Se observa que hay problemas de multicolinealidad, homogenidad de la varianza, residuales no normales y linealidad.
Y en el grafico de bandas simuladas se alcanzan a salir varios puntos, esto debido a los problemas mencionados anteriormente.
d)
Interprete las estimaciones de los parámetros excepto el intercepto.
La media de cada individuo, segun el modelo es de la forma, donde la referencia base es es ser hombre:
\[ \mu_i=\beta_0+\beta_1*Sexo_i+\beta_2*Edad_i+\beta_3*(Sexo_i*Edad_i) \] Entonces, un aumento de una unidad en la edad de gestación para un hombre, aumenta en 111.98 unidades el valor promedio del peso.
Y un aumento en una unidad en la edad de gestación para una mujer aumenta en 111.98+18.42=130.4 unidades el valor promedio del peso.
e)
Estime el peso esperado de una bebita con 40 semanas de gestación.
bebe = data.frame(Sexo="Mujer",Edad=40)
predict(fit1,newdata =bebe )## 1
## 3074.333Punto 2
Primero se cargan los datos.
datos2 <- read.table(header = TRUE, text = "
Escala Senilidad
9 1
13 1
6 1
8 1
10 1
4 1
14 1
8 1
11 1
7 1
9 1
7 1
5 1
14 1
13 0
16 0
10 0
12 0
11 0
14 0
15 0
18 0
7 0
16 0
9 0
9 0
11 0
13 0
15 0
13 0
10 0
11 0
6 0
17 0
14 0
19 0
9 0
11 0
14 0
10 0
16 0
16 0
14 0
13 0
13 0
9 0
15 0
10 0
11 0
12 0
4 0
14 0
20 0
10 0
")
datos2$Senilidad <- as.factor(datos2$Senilidad)a)
Ajuste modelos de respuesta binaria con funciones de enlace logit, probit y complemento log-log. Use el desvio y el criterio AIC para seleccionar el “mejor” modelo.
modelo_logit <- glm(Senilidad ~ Escala,
data = datos2,
family = binomial(link = "logit"))
summary(modelo_logit)##
## Call:
## glm(formula = Senilidad ~ Escala, family = binomial(link = "logit"),
## data = datos2)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.4040 1.1918 2.017 0.04369 *
## Escala -0.3235 0.1140 -2.838 0.00453 **
## ---
## 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: 51.017 on 52 degrees of freedom
## AIC: 55.017
##
## Number of Fisher Scoring iterations: 5modelo_probit <- glm(Senilidad ~ Escala,
data = datos2,
family = binomial(link = "probit"))
summary(modelo_probit)##
## Call:
## glm(formula = Senilidad ~ Escala, family = binomial(link = "probit"),
## data = datos2)
##
## 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: 5modelo_cloglog <- glm(Senilidad ~ Escala,
data = datos2,
family = binomial(link = "cloglog"))
summary(modelo_cloglog)##
## Call:
## glm(formula = Senilidad ~ Escala, family = binomial(link = "cloglog"),
## data = datos2)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.42612 0.81686 1.746 0.0808 .
## Escala -0.25075 0.08389 -2.989 0.0028 **
## ---
## 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: 51.311 on 52 degrees of freedom
## AIC: 55.311
##
## Number of Fisher Scoring iterations: 6El mejor modelo tanto por AIC, como por devianza es el modelo con el enlace probit.
b)
Use las estadísticas de Wald y razón de verosimilitudes para evaluar, al nivel de significancia aproximado de 5%, si la probabilidad de presentar síntomas de senilidad depende del puntaje en la escala de inteligencia de Wechsler.
modelo_nulo <- glm(Senilidad ~ 1,
data = datos2,
family = binomial(link = "probit"))
anova(modelo_nulo, modelo_probit, test = "LRT")## Analysis of Deviance Table
##
## Model 1: Senilidad ~ 1
## Model 2: Senilidad ~ Escala
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 53 61.806
## 2 52 50.984 1 10.823 0.001003 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En el summary del punto anterior se observa que el test de Wald con una significancia del 5%, rechaza la hipotesis nula, es decir el parametro es estadiscamente significativo.
Y para el test de razón de verosimilitud se concluye que el modelo con la variable es diferente al modelo nulo.
c)
Realice el análisis de diagnóstico. Comente.
check_model(modelo_probit)
set.seed(123)
envelope(modelo_probit, rep=1000)
No se observa nada critico en el diagnostico del modelo, mas allá de una cola pesada en el qqplot de check_model, pero con las bandas simuladas no se parecia este problema, porque caen dentro de las bandas.
d)
Interprete las estimaciones de los parámetros excepto el intercepto.
Como la función de enlace que se esta usando es la probit, dificulta demasiado la interpretación de los paramaetros, lo uncio que se puede asegurar es que como \(\hat\beta_1\) es negativo a medida que aumenta la escala de inteligencia de Wechsler (que solo toma valores positivos), disminuye la probabilidad promedio de tener senilidad.
e)
Cuál es la probabilidad estimada que una persona con 15 puntos en la escala de Wechsler tenga síntomas de senilidad?
persona_new = data.frame(Escala=15)
predict(modelo_probit,newdata = persona_new, type ="response" )## 1
## 0.07577691Punto 3
Se cargan los datos y como se ajusta un modelo binomial, se crea el vector “respuesta”.
datos3 <- read.table(header = TRUE, text = "
Dosis Expuestos Muertos
0.41 50 6
0.58 48 16
0.71 46 24
0.89 49 42
1.01 50 44
")
respuesta = cbind(datos3$Muertos, datos3$Expuestos-datos3$Muertos)
dosis = datos3$Dosisa)
Realice un gráfico de la tasa de mortalidad de los insectos como funciónn de la concentración del plaguicida. Comente.
plot(y=datos3$Muertos/datos3$Expuestos,x=datos3$Dosis, ylab="Tasa de Mortalidad", xlab="Dosis" )
Se aprecia que a medida que se aumenta la dosis, la tasa de mortalidad aumenta y aparenta ser una relación lineal.
b)
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(datos3$Muertos, datos3$Expuestos-datos3$Muertos)
modelo_logit <- glm(respuesta ~ dosis, family=binomial)
modelo_probit <- glm(respuesta ~ dosis, family=binomial(link="probit"))
modelo_cloglog <- glm(respuesta ~ dosis, family=binomial(link="cloglog"))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
El modelo que mejor predice la tasa, osea el que mas se acerca a los puntos negros es la linea azul, el logit.
c)
Compare los tres modelos ajustados usando los valores del desvio y del AIC. Seleccione el “mejor” modelo.
respuesta = cbind(datos3$Muertos, datos3$Expuestos-datos3$Muertos)
modelo_logit <- glm(respuesta ~ dosis, family=binomial)
summary(modelo_logit)##
## Call:
## glm(formula = respuesta ~ dosis, family = binomial)
##
## 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: 4modelo_probit <- glm(respuesta ~ dosis, family=binomial(link="probit"))
summary(modelo_probit)##
## Call:
## glm(formula = respuesta ~ dosis, family = binomial(link = "probit"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.8594 0.3481 -8.214 <2e-16 ***
## dosis 4.1691 0.4752 8.774 <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: 96.6881 on 4 degrees of freedom
## Residual deviance: 1.6386 on 3 degrees of freedom
## AIC: 24.859
##
## Number of Fisher Scoring iterations: 4modelo_cloglog <- glm(respuesta ~ dosis, family=binomial(link="cloglog"))
summary(modelo_cloglog)##
## Call:
## glm(formula = respuesta ~ dosis, family = binomial(link = "cloglog"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.5119 0.4304 -8.159 3.38e-16 ***
## dosis 4.4259 0.5244 8.441 < 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: 96.6881 on 4 degrees of freedom
## Residual deviance: 4.0355 on 3 degrees of freedom
## AIC: 27.256
##
## Number of Fisher Scoring iterations: 4El modelo con mejores métricas es el logit.
d)
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.
summary(modelo_logit)##
## Call:
## glm(formula = respuesta ~ dosis, family = binomial)
##
## 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: 4modelo_nulo <- glm(respuesta ~ 1, family=binomial)
anova(modelo_logit,modelo_nulo)## Analysis of Deviance Table
##
## Model 1: respuesta ~ dosis
## Model 2: respuesta ~ 1
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 3 1.348
## 2 4 96.688 -1 -95.34 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En el test de wald el parametro es estadísticamente significativo y en el de razón de verosimilitudes es difernete respecto al modelo nulo.
e)
Realice el análisis de diagnóstico. Comente.
check_model(modelo_probit)## Some of the variables were in matrix-format - probably you used
## `scale()` on your data?
## If so, and you get an error, please try `datawizard::standardize()` to
## standardize your data.
## Some of the variables were in matrix-format - probably you used
## `scale()` on your data?
## If so, and you get an error, please try `datawizard::standardize()` to
## standardize your data.
set.seed(123)
envelope(modelo_probit, rep=1000)
Es un poco dificil evaluar supuestos con pocos datos, pero en general no se detecta ningún problema.
f)
Estime la dosis letal 50 % y calcule para ella un intervalo de confianza de aproximadamente 95 %.
Recordemos que la dosis letal debe cumplir que:
\[ log(\frac{0.5}{1-0.5})=\beta_0+\beta_1*Dosis_{Let} \] \[ 0=\beta_0+\beta_1*Dosis_{Let} \] \[ Dosis_{Let}=\frac{-\beta_0}{\beta_1} \] Y los intervalos se calculan por medio de una función de R porque la dosis letal no posee una distribución trivial.
-4.8387/ -7.0681## [1] 0.6845828dose.p(modelo_logit)## Dose SE
## p = 0.5: 0.6845816 0.02271034c(0.6845816 -1.96*0.02271034,0.6845816 +1.96*0.02271034)## [1] 0.6400693 0.7290939Punto 4
Se cargan los datos.
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)
)a)
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.
###Modelos higado
dosisH =datos_higado$Dosis
respuestaH = cbind(datos_higado$Enfermos, datos_higado$Expuestos-datos_higado$Enfermos)
modelo_logitH <- glm(respuestaH ~ dosisH, family=binomial)
summary(modelo_logitH)##
## Call:
## glm(formula = respuestaH ~ dosisH, family = binomial)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.5083 0.1478 -30.501 <2e-16 ***
## dosisH 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: 4modelo_poisH = glm(Enfermos ~ Dosis+ offset(log(Expuestos)), family=poisson(link="log"),
data=datos_higado)
summary(modelo_poisH)##
## 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#### Modelos vejiga
dosisV =datos_vejiga$Dosis
respuestaV = cbind(datos_vejiga$Enfermos, datos_vejiga$Expuestos-datos_vejiga$Enfermos)
modelo_logitV <- glm(respuestaV ~ dosisV, family=binomial)
summary(modelo_logitV)##
## Call:
## glm(formula = respuestaV ~ dosisV, family = binomial)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.4324 0.5831 -12.75 <2e-16 ***
## dosisV 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: 5modelo_poisV = glm(Enfermos ~ Dosis+ offset(log(Expuestos)), family=poisson(link="log"),
data=datos_vejiga)
summary(modelo_poisV)##
## 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: 5El modelo seleccionado en ambos casos es el logit, ya que presenta un mejor AIC. En este contexto, los modelos son comparables porque el modelo Poisson incluye un offset, lo que permite interpretarlo en términos de tasas o proporciones, haciéndolo coherente con la comparación frente al binomial. Además, para este caso clínico resulta útil estimar cantidades como la dosis letal.
b)
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´ogeno afecta la probabilidad/incidencia de cáncer en los ratones.
modelo_logitVnu <- glm(respuestaV ~ 1, family=binomial)
modelo_logitHnu <- glm(respuestaH ~ 1, family=binomial)
anova(modelo_logitVnu,modelo_logitV)## Analysis of Deviance Table
##
## Model 1: respuestaV ~ 1
## Model 2: respuestaV ~ dosisV
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 7 214.94
## 2 6 11.45 1 203.5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(modelo_logitHnu,modelo_logitH)## Analysis of Deviance Table
##
## Model 1: respuestaH ~ 1
## Model 2: respuestaH ~ dosisH
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 7 77.603
## 2 6 4.522 1 73.081 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En los dos casos, el parametro es estadisticamente significativo y los modelos son diferentes a su respectivo modelo nulo.
e)
Realice el análisis de diagnóstico. Comente.
check_model(modelo_logitH)## Some of the variables were in matrix-format - probably you used
## `scale()` on your data?
## If so, and you get an error, please try `datawizard::standardize()` to
## standardize your data.
## Some of the variables were in matrix-format - probably you used
## `scale()` on your data?
## If so, and you get an error, please try `datawizard::standardize()` to
## standardize your data.
set.seed(123)
envelope(modelo_logitH, rep = 10000)
check_model(modelo_logitV)## Some of the variables were in matrix-format - probably you used
## `scale()` on your data?
## If so, and you get an error, please try `datawizard::standardize()` to
## standardize your data.
## Some of the variables were in matrix-format - probably you used
## `scale()` on your data?
## If so, and you get an error, please try `datawizard::standardize()` to
## standardize your data.
set.seed(123)
envelope(modelo_logitV, rep = 10000)
Hay algunos individuos con alto leverage y residuos estandarizados, es decir pueden ser observaciones influyentes. Del resto todo se observa bien.
#Hay unos valores influyentes ver si cambian las estimaciones
source("C:\\Users\\User\\Desktop\\MLG\\Diapositivas\\GLM_Count\\macros.txt")
CDH <- cooks.distance(modelo_logitH)
plot(CDH,type="h")
text(which(CDH>1), CDH[which(CDH>1)], which(CDH>1))
case.deletion_glm(modelo_logitH, subset = -c(8))##
## Estimate Std. Error Pr(>|z|) | Estimate* Std. Error* Pr(>|z|)*
## (Intercept) -4.5083 0.1478 0 | -4.7436 0.1987 0
## dosisH 1.7624 0.1864 0 | 2.3175 0.3486 0
## | Change(%)
## (Intercept) | -5.219
## dosisH | 31.495
##
## (*) estimates, standard errors and p-values obtained using the specified subset of individuals.CDV <- cooks.distance(modelo_logitV)
plot(CDV,type="h")
text(which(CDV>1), CDV[which(CDV>1)], which(CDV>1))
case.deletion_glm(modelo_logitV, subset = -c(2,7))##
## Estimate Std. Error Pr(>|z|) | Estimate* Std. Error* Pr(>|z|)*
## (Intercept) -7.4324 0.5831 0 | -7.7809 1.0604 0
## dosisV 7.8752 0.7821 0 | 8.0197 1.5581 0
## | Change(%)
## (Intercept) | -4.689
## dosisV | 1.835
##
## (*) estimates, standard errors and p-values obtained using the specified subset of individuals.case.deletion_glm(modelo_logitV, subset = -c(7))##
## Estimate Std. Error Pr(>|z|) | Estimate* Std. Error* Pr(>|z|)*
## (Intercept) -7.4324 0.5831 0 | -6.9512 0.6311 0
## dosisV 7.8752 0.7821 0 | 6.9038 1.0016 0
## | Change(%)
## (Intercept) | 6.474
## dosisV | -12.335
##
## (*) estimates, standard errors and p-values obtained using the specified subset of individuals.case.deletion_glm(modelo_logitV, subset = -c(2))##
## Estimate Std. Error Pr(>|z|) | Estimate* Std. Error* Pr(>|z|)*
## (Intercept) -7.4324 0.5831 0 | -8.3187 0.8622 0
## dosisV 7.8752 0.7821 0 | 8.9208 1.0840 0
## | Change(%)
## (Intercept) | -11.92
## dosisV | 13.28
##
## (*) estimates, standard errors and p-values obtained using the specified subset of individuals.Se aprecia que quitando los individuos no hay cambio de signo y los cambios de magnitud de los parametros no son muy grandes
d)
Interprete las estimaciones de los par´ametros excepto los interceptos.
Los aumentos porcentuales son de la forma
\[ exp(\beta_j)-1 \]
#Higado
exp(0.1* 1.7624)-1## [1] 0.1927243#Vejiga
exp(0.1*7.8752)-1## [1] 1.197939Un aumento de 0.1 unidades en las dosis, aumenta los chances promedio de tener cancer de higado y vejiga en 19% y 119% respectivamente.
Punto 5
Se cargan los datos
datos5 <- data.frame(
Tiempo_Modo_1 = c(33.3, 52.2, 64.7, 137.0, 125.9, 116.3, 131.7, 85.0, 91.9),
Tiempo_Modo_2 = c(25.3, 14.4, 32.5, 20.5, 97.6, 53.6, 56.6, 87.3, 47.8),
Num_Fallas = c(15, 9, 14, 24, 27, 27, 23, 18, 22)
)a)
Ajuste modelos de respuesta poisson con funci´on de enlace logarítmica, identidad y raiz 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.
fit_log = glm(Num_Fallas ~ Tiempo_Modo_1+Tiempo_Modo_2, family=poisson(link="log"),
data=datos5)
fit_iden = glm(Num_Fallas ~ Tiempo_Modo_1+Tiempo_Modo_2, family=poisson(link="identity"),
data=datos5)
fit_sqrt = glm(Num_Fallas ~ Tiempo_Modo_1+Tiempo_Modo_2, family=poisson(link="sqrt"),
data=datos5)
summary(fit_log)##
## Call:
## glm(formula = Num_Fallas ~ Tiempo_Modo_1 + Tiempo_Modo_2, family = poisson(link = "log"),
## data = datos5)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.175168 0.255456 8.515 < 2e-16 ***
## Tiempo_Modo_1 0.007015 0.002429 2.888 0.00387 **
## Tiempo_Modo_2 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: 4summary(fit_iden)##
## Call:
## glm(formula = Num_Fallas ~ Tiempo_Modo_1 + Tiempo_Modo_2, family = poisson(link = "identity"),
## data = datos5)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.99773 3.63545 1.650 0.09899 .
## Tiempo_Modo_1 0.12081 0.04578 2.639 0.00832 **
## Tiempo_Modo_2 0.05459 0.06356 0.859 0.39037
## ---
## 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.1971 on 6 degrees of freedom
## AIC: 53.254
##
## Number of Fisher Scoring iterations: 6summary(fit_sqrt)##
## Call:
## glm(formula = Num_Fallas ~ Tiempo_Modo_1 + Tiempo_Modo_2, family = poisson(link = "sqrt"),
## data = datos5)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.761295 0.485366 5.689 1.28e-08 ***
## Tiempo_Modo_1 0.014771 0.005245 2.816 0.00486 **
## Tiempo_Modo_2 0.005820 0.006743 0.863 0.38807
## ---
## 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.0341 on 6 degrees of freedom
## AIC: 53.091
##
## Number of Fisher Scoring iterations: 4El mejor modelo por AIC, es el de enlace log.
b y c)
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´pnicas depende del tiempo de funcionamiento en el modo 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´pnicas depende del tiempo de funcionamiento en el modo 2.
summary(fit_log)##
## Call:
## glm(formula = Num_Fallas ~ Tiempo_Modo_1 + Tiempo_Modo_2, family = poisson(link = "log"),
## data = datos5)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.175168 0.255456 8.515 < 2e-16 ***
## Tiempo_Modo_1 0.007015 0.002429 2.888 0.00387 **
## Tiempo_Modo_2 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: 4fit_mod2 = glm(Num_Fallas ~ Tiempo_Modo_2, family=poisson(link="log"),
data=datos5)
anova(fit_log, fit_mod2)## Analysis of Deviance Table
##
## Model 1: Num_Fallas ~ Tiempo_Modo_1 + Tiempo_Modo_2
## Model 2: Num_Fallas ~ Tiempo_Modo_2
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 6 4.0033
## 2 7 12.3561 -1 -8.3528 0.003851 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit_mod1 = glm(Num_Fallas ~ Tiempo_Modo_1, family=poisson(link="log"),
data=datos5)
anova(fit_log, fit_mod1)## Analysis of Deviance Table
##
## Model 1: Num_Fallas ~ Tiempo_Modo_1 + Tiempo_Modo_2
## Model 2: Num_Fallas ~ Tiempo_Modo_1
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 6 4.0033
## 2 7 4.8078 -1 -0.80451 0.3697El parametro del modo 2, no es significativo y ademas la prueba de razón de verosimilitud, dice que el modelo con el modo 1 y 2, es estadisitacamente igual al que solo tiene modo 1. Por lo tanto se elimina esa variable del modelo.
d)
Realice el análisis de diagnóstico. Comente.
check_model(fit_mod1)
set.seed(123)
envelope(fit_mod1, rep = 10000)
Se observa tanto en el check model y el envelope, que no hay mal especificación del modelo como muchos 0 o sobredispersión o subdispersión.
e)
Interprete las estimaciones de los parámetros excepto el intercepto.
summary(fit_mod1)##
## Call:
## glm(formula = Num_Fallas ~ Tiempo_Modo_1, family = poisson(link = "log"),
## data = datos5)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.237196 0.243053 9.205 < 2e-16 ***
## Tiempo_Modo_1 0.007705 0.002264 3.403 0.000667 ***
## ---
## 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.8078 on 7 degrees of freedom
## AIC: 51.865
##
## Number of Fisher Scoring iterations: 4La interpretación se hace de la siguiente manera:
\[ exp(\beta_j)-1 \]
exp(0.007705)-1## [1] 0.00773476Por cada aumento de una unidad en el tiempo de funcionamiento del modo 1, el numero de fallas promedio aumenta en 0.77%.
Punto 6
Cargamos los datos.
datos6 <- data.frame(
Tiempo_Exposicion = c(0, 15, 30, 45, 60),
Num_Lesiones = c(271, 108, 59, 29, 12)
)a)
Realice un gráfico del número de lesiones como funci´on del tiempo de exposición. Comente.
plot(datos6$Tiempo_Exposicion,datos6$Num_Lesiones,xlab="Tiempo de Expo",ylab="# Lesiones")
Se puede apreciar que a mayor tiempo de exposición hay un menor numero de lesiones, y no aparenta ser una relación lineal.
b)
Ajuste modelos de respuesta poisson con función de enlace logarítmica, identidad y raiz cuadrada. Use el desvio y el criterio AIC para seleccionar el “mejor” modelo.
fit_log = glm(Num_Lesiones ~ Tiempo_Exposicion, family=poisson(link="log"),
data=datos6)
fit_iden = glm(Num_Lesiones ~ Tiempo_Exposicion, family=poisson(link="identity"),
data=datos6)
fit_sqrt = glm(Num_Lesiones ~ Tiempo_Exposicion, family=poisson(link="sqrt"),
data=datos6)
summary(fit_log)##
## Call:
## glm(formula = Num_Lesiones ~ Tiempo_Exposicion, family = poisson(link = "log"),
## data = datos6)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.571340 0.056714 98.24 <2e-16 ***
## Tiempo_Exposicion -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: 3summary(fit_iden)##
## Call:
## glm(formula = Num_Lesiones ~ Tiempo_Exposicion, family = poisson(link = "identity"),
## data = datos6)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 183.9983 8.8115 20.88 <2e-16 ***
## Tiempo_Exposicion -2.9399 0.1606 -18.31 <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.131 on 4 degrees of freedom
## Residual deviance: 74.184 on 3 degrees of freedom
## AIC: 107.61
##
## Number of Fisher Scoring iterations: 13summary(fit_sqrt)##
## Call:
## glm(formula = Num_Lesiones ~ Tiempo_Exposicion, family = poisson(link = "sqrt"),
## data = datos6)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 14.94189 0.38730 38.58 <2e-16 ***
## Tiempo_Exposicion -0.20620 0.01054 -19.56 <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.131 on 4 degrees of freedom
## Residual deviance: 25.947 on 3 degrees of freedom
## AIC: 59.375
##
## Number of Fisher Scoring iterations: 7El modelo con mejor AIC, es el que tiene enlace log.
c)
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(fit_log)##
## Call:
## glm(formula = Num_Lesiones ~ Tiempo_Exposicion, family = poisson(link = "log"),
## data = datos6)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.571340 0.056714 98.24 <2e-16 ***
## Tiempo_Exposicion -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: 3fit_nu = glm(Num_Lesiones ~ 1, family=poisson(link="log"),
data=datos6)
anova(fit_log,fit_nu)## Analysis of Deviance Table
##
## Model 1: Num_Lesiones ~ Tiempo_Exposicion
## Model 2: Num_Lesiones ~ 1
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 3 2.29
## 2 4 413.13 -1 -410.84 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El parametro es signifcativo segun el test de Wald y el test de razón de verosimilitud dice que hay diferencias entre el modelo nulo y el que tiene la variable explicativa.
d)
Realice el análisis de diagnóstico. Comente.
check_model(fit_log)
set.seed(123)
envelope(fit_log, rep = 10000)
Hay mus pocos datos, por lo que el diagnostico es muy poco preciso, pero se observan falencias en la homogenidad de la varianza y la mal especifiación de muchos 0 o sobredispersión. Pese a esto el grafico de bandas simuladas se ve bien
e)
Interprete las estimaciones de los parámetros excepto el intercepto.
La interpretación se hace de la siguiente manera:
\[ exp(\beta_j)-1 \]
exp(-0.051326)-1## [1] -0.05003107Por cada aumento de una unidad en el tiempo de exposición, el numero de lesiones promedio disminuye en 5%.
f)
Estime el número esperado de lesiones causadas por el virus después de su exposición a rayos X durante 30 minutos.
new=data.frame(Tiempo_Exposicion=30)
predict(fit_log,newdata = new, type="response")## 1
## 56.34955Punto 7
Se cargan los datos.
datos7 <- data.frame(
Ciudad = c(rep("St_Paul", 7), rep("Forth_Worth", 7)),
Edad = c("15-24", "25-34", "35-44", "45-54", "55-64", "65-74", "75+",
"15-24", "25-34", "35-44", "45-54", "55-64", "65-74", "75+"),
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)
)
# Mostrar la estructura para verificar
datos7$Ciudad = as.factor(datos7$Ciudad)
datos7$Edad = as.factor(datos7$Edad)a)
Ajuste un modelo de respuesta Poisson con función de enlace logarítmica y las covariables cualitativas ciudad y edad. Comente.
fit_log = glm(Casos ~ Ciudad+Edad+offset(log(Poblacion)), family=poisson(link="log"),
data=datos7)
summary(fit_log)##
## Call:
## glm(formula = Casos ~ Ciudad + Edad + offset(log(Poblacion)),
## family = poisson(link = "log"), data = datos7)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.83938 0.44754 -24.220 < 2e-16 ***
## CiudadSt_Paul -0.85269 0.05962 -14.302 < 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 ***
## ---
## 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: 4Se aprecia que la devianza residual y los grados de libertad, son cercanos lo que indica un buen ajuste y ademas todos los parametros son significativos.
b)
fit_nu = glm(Casos ~ Ciudad+offset(log(Poblacion)), family=poisson(link="log"),
data=datos7)
anova(fit_log,fit_nu)## Analysis of Deviance Table
##
## Model 1: Casos ~ Ciudad + Edad + offset(log(Poblacion))
## Model 2: Casos ~ Ciudad + offset(log(Poblacion))
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 6 5.21
## 2 12 1957.76 -6 -1952.6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(fit_log)##
## Call:
## glm(formula = Casos ~ Ciudad + Edad + offset(log(Poblacion)),
## family = poisson(link = "log"), data = datos7)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.83938 0.44754 -24.220 < 2e-16 ***
## CiudadSt_Paul -0.85269 0.05962 -14.302 < 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 ***
## ---
## 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: 4El test de Wald dice que el parametro de la ciudad es significativo, es decir cancer de piel depende de la ciudad y el test de razón de verosimilitud dice que el modelo solo con edad es significativamente diferente que el que tiene las dos variables.
c)
Realice el análisis de diagnóstico. Comente.
check_model(fit_log)
set.seed(123)
envelope(fit_log, rep = 10000)
Al parecer hay mal espicifación del modelo, por el grafico de los residuales parece sobredispersión.
d)
La incidencia de cancer de piel es mayor en Fort Worth que en St Paul? Justifique.
Como el coeficiente de St.Paul es negativo, indica que a comparación de Forth Worth la tasas de cáncer de piel en pormedio es menor ahí.
e)
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 10.000 personas.
new=data.frame(Ciudad="St_Paul",Edad="75+",Poblacion=10000)
predict(fit_log,newdata = new,type="response")## 1
## 40.20308Punto 8
Cargamos los datos
data(wafer)a)
Ajuste modelos de respuesta Gama con funciones de enlace identidad, logarítmica y recíproca, y los cuatros factores como covariables. Use el criterio AIC para seleccionar el “mejor” modelo.
fit_iden = glm(resist ~ x1+x2+x3+x4, family=Gamma(link="identity"),
data=wafer)
fit_log = glm(resist ~ x1+x2+x3+x4, family=Gamma(link="log"),
data=wafer)
fit_inv = glm(resist ~ x1+x2+x3+x4, family=Gamma(link="inverse"),
data=wafer)
summary(fit_iden)##
## Call:
## glm(formula = resist ~ x1 + x2 + x3 + x4, family = Gamma(link = "identity"),
## data = wafer)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 235.43 14.21 16.572 3.97e-09 ***
## x1+ 27.83 12.24 2.274 0.044017 *
## x2+ -65.74 12.68 -5.184 0.000302 ***
## x3+ 36.94 12.32 2.999 0.012102 *
## x4+ -11.88 12.15 -0.978 0.349294
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.0124071)
##
## Null deviance: 0.69784 on 15 degrees of freedom
## Residual deviance: 0.14009 on 11 degrees of freedom
## AIC: 154.84
##
## Number of Fisher Scoring iterations: 6summary(fit_log)##
## Call:
## glm(formula = resist ~ x1 + x2 + x3 + x4, family = Gamma(link = "log"),
## data = wafer)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.44552 0.05856 92.983 < 2e-16 ***
## x1+ 0.12115 0.05238 2.313 0.041090 *
## x2+ -0.30049 0.05238 -5.736 0.000131 ***
## x3+ 0.17979 0.05238 3.432 0.005601 **
## x4+ -0.05757 0.05238 -1.099 0.295248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.01097542)
##
## Null deviance: 0.69784 on 15 degrees of freedom
## Residual deviance: 0.12418 on 11 degrees of freedom
## AIC: 152.91
##
## Number of Fisher Scoring iterations: 4summary(fit_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: 4El mejor modelo por AIC es el inverso, pero hay varias variables que no son significativas.
b)
Evalúe, al nivel de significancia aproximado de 5 %, cual(es) factor(es) afecta(n) la resistencia de las placas de los semiconductores.
Las variables significativas son x2 y x3, entonces se ajusta un modelo con esas variables.
fit_inv3 = glm(resist ~ x2+x3, family=Gamma(link="inverse"),
data=wafer)
summary(fit_inv3)##
## Call:
## glm(formula = resist ~ x2 + x3, family = Gamma(link = "inverse"),
## data = wafer)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0042377 0.0002180 19.437 5.44e-11 ***
## x2+ 0.0013498 0.0002613 5.166 0.000182 ***
## x3+ -0.0008183 0.0002540 -3.222 0.006683 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.01348383)
##
## Null deviance: 0.69784 on 15 degrees of freedom
## Residual deviance: 0.17980 on 13 degrees of freedom
## AIC: 154.84
##
## Number of Fisher Scoring iterations: 4c)
Realice el análisis de diagnóstico. Comente.
check_model(fit_inv3)
set.seed(123)
envelope(fit_inv3, rep = 10000)
No se aprecia ningun indicio de mal ajuste.