Taller 3_Modelos Lineales Generalizados
JOHN JAIRO PRADO
2024-10-26
1. En la base de datos blocks, se dan datos de un experimento donde se pidió a los niños que construyeran torres con bloques cúbicos y cilíndricos como tan alto como pudieron, el número de bloques utilizados y el tiempo empleado es la variable de interés.
(1a) Ajuste el “mejor” modelo de regresión con el tiempo que tardaron los niños como variable de interés. Muestre las medidas de bondad de ajuste.
Descripcion del Ejercicio
Child: the child; an identifier from A to Y. El niño; un identificador de A a Y la cantidad de bloques que el niño podría apilar con éxito; un vector numérico
Number: the number of blocks the child could successfully stack; a numeric vector
Time: the time in seconds taken for the children to make their stack of blocks; a numeric vector. El tiempo en segundos que les tomó a los niños hacer su pila de bloques; un vector numérico.
Trial: the trial number on which the data were gathered (see Details); a factor with levels 1 and 2. El número de ensayo sobre el cual se recopilaron los datos (ver Detalles); un factor con niveles 1 y 2
Shape: the shape of the blocks being stacked; a factor with levels Cube and Cylinder. la forma de los bloques que se apilan; un factor con niveles Cubo y Cilindro
Age: the age of the child in completed years; a numeric vector. La edad del niño en años cumplidos; un vector numérico
suppressPackageStartupMessages(library(AER))
suppressPackageStartupMessages(library(GLMsData))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(MASS, exclude = "select"))
suppressPackageStartupMessages(library(plotly))data(blocks)
glimpse(blocks)Rows: 100
Columns: 6
$ Child <fct> A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, T, U, …
$ Number <int> 11, 9, 8, 9, 10, 13, 10, 7, 6, 12, 6, 9, 8, 7, 9, 9, 8, 7, 8, 4…
$ Time <dbl> 30.0, 19.0, 18.6, 23.0, 29.0, 178.0, 42.0, 20.2, 12.0, 39.2, 14…
$ Trial <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ Shape <fct> Cube, Cube, Cube, Cube, Cube, Cube, Cube, Cube, Cube, Cube, Cub…
$ Age <dbl> 4.67, 5.00, 4.42, 4.33, 4.33, 4.83, 4.42, 5.00, 3.58, 4.00, 4.5…hist(blocks$Time)
Se observa una asimetría positiva.
blocks$Child_dummy <- ifelse(blocks$Child %in% LETTERS[1:12], 1, 0)
# Crear la variable dummy para 'Number'
blocks$Number_dummy <- ifelse(blocks$Number >= 3 & blocks$Number <= 7, 1, 0)
blocks$Trial_dummy=ifelse(blocks$Trial==2,1,0)
blocks$Shape_dummy=ifelse(blocks$Shape=="Cube",1,0)
blocks$Age_dummy <- ifelse(blocks$Age >= 2.42 & blocks$Age <= 4.00, 1, 0)
# Eliminar columnas asignando a NULL
blocks$Child <- NULL
blocks$Number <- NULL
blocks$Trial <- NULL
blocks$Shape <- NULL
blocks$Age <- NULL
# Renombrar columnas usando colnames()
colnames(blocks)[colnames(blocks) == "Child_dummy"] <- "Child"
colnames(blocks)[colnames(blocks) == "Number_dummy"] <- "Number"
colnames(blocks)[colnames(blocks) == "Trial_dummy"] <- "Trial"
colnames(blocks)[colnames(blocks) == "Shape_dummy"] <- "Shape"
colnames(blocks)[colnames(blocks) == "Age_dummy"] <- "Age"
head(blocks) Time Child Number Trial Shape Age
1 30.0 1 0 0 1 0
2 19.0 1 0 0 1 0
3 18.6 1 0 0 1 0
4 23.0 1 0 0 1 0
5 29.0 1 0 0 1 0
6 178.0 1 0 0 1 0Verificar si es Gamma o Inversa Gaussiana: Si el coeficiente de variacion e mayor a 0.4 la inversa gaussiana es una buena opcion de modelado. y si queda entre esos rangos hay que realizar un test de hipotesis para verificar
modelo_ga_blocks=glm(Time~.,data=blocks,family=Gamma(link="log"))
modelo_inv_blocks=glm(Time~.,data=blocks,family=inverse.gaussian(link="log"))
log_likelihood = logLik(modelo_ga_blocks)/logLik(modelo_inv_blocks); log_likelihood'log Lik.' 1.013 (df=7)Es importante colocar (link=“log”) para que los resultados esten en esa funcion
Si log_likelihood > 0 es inversamente gaussiana Si log_likelihood < 0 es gamma
CONCLUSION: Como El logaritmo del cociente de las verosimilitudes (log_likelihood) es 1.013 se debe utilizar el modelo de regresion inversamente Gaussiana
Se comprueba con el test de Hipotesis:
modelo=lm(Time~.,data=blocks)
shapiro.test(modelo$residuals)
Shapiro-Wilk normality test
data: modelo$residuals
W = 0.69, p-value = 4e-13Hipótesis nula (H0): Los datos (en este caso, los residuos) provienen de una distribución normal. Hipótesis alternativa (H1): Los datos no provienen de una distribución normal.
Como el p-valor es menor que un nivel de significancia predefinido (como 0.05), se rechaza la hipótesis nula.
Seleccion del Mejor Subconjunto de Variables
step(modelo_inv_blocks)Start: AIC=801.9
Time ~ Child + Number + Trial + Shape + Age
Df Deviance AIC
- Trial 1 1.28 800
- Child 1 1.29 800
- Shape 1 1.31 802
<none> 1.28 802
- Age 1 1.36 804
- Number 1 1.75 828
Step: AIC=800.1
Time ~ Child + Number + Shape + Age
Df Deviance AIC
- Child 1 1.29 798
- Shape 1 1.31 800
<none> 1.28 800
- Age 1 1.36 803
- Number 1 1.76 826
Step: AIC=798.6
Time ~ Number + Shape + Age
Df Deviance AIC
- Shape 1 1.32 798
<none> 1.29 799
- Age 1 1.39 802
- Number 1 1.81 828
Step: AIC=798.4
Time ~ Number + Age
Df Deviance AIC
<none> 1.32 798
- Age 1 1.42 802
- Number 1 2.03 836
Call: glm(formula = Time ~ Number + Age, family = inverse.gaussian(link = "log"),
data = blocks)
Coefficients:
(Intercept) Number Age
3.746 -0.970 0.342
Degrees of Freedom: 99 Total (i.e. Null); 97 Residual
Null Deviance: 2.03
Residual Deviance: 1.32 AIC: 798mejor_modelo_inv_blocks = glm(Time ~ Number + Age, data = blocks, family=inverse.gaussian(link="log"))
summary(mejor_modelo_inv_blocks)
Call:
glm(formula = Time ~ Number + Age, family = inverse.gaussian(link = "log"),
data = blocks)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.746 0.163 22.96 < 2e-16 ***
Number -0.970 0.175 -5.54 2.6e-07 ***
Age 0.342 0.134 2.55 0.012 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for inverse.gaussian family taken to be 0.01808)
Null deviance: 2.0308 on 99 degrees of freedom
Residual deviance: 1.3156 on 97 degrees of freedom
AIC: 798.4
Number of Fisher Scoring iterations: 6Medidas de Bondad de Ajuste
Pseudo R^2
library(DescTools)
Adjuntando el paquete: 'DescTools'The following object is masked from 'package:car':
RecodePseudoR2(modelo_inv_blocks, 'Nagelkerke')Nagelkerke
0.3678 Nota: x < 0.2 El estadistico Nagelkerke recomienda una mala explicacion del modelo de la variable de interes.
0.21 > x < 0.4 El estadistico Nagelkerke recomienda que es una buena explicacion del modelo de la variable de interes.
x > 0.41 El estadistico Nagelkerke recomienda una excelente explicacion del modelo de la variable de interes.
El modelo de regresion inversamente gaussiana explica de manera buena el total de tiempo en segundos que les tomó a los niños hacer su pila de bloques
Prueba de la Deviance
summary(mejor_modelo_inv_blocks)
Call:
glm(formula = Time ~ Number + Age, family = inverse.gaussian(link = "log"),
data = blocks)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.746 0.163 22.96 < 2e-16 ***
Number -0.970 0.175 -5.54 2.6e-07 ***
Age 0.342 0.134 2.55 0.012 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for inverse.gaussian family taken to be 0.01808)
Null deviance: 2.0308 on 99 degrees of freedom
Residual deviance: 1.3156 on 97 degrees of freedom
AIC: 798.4
Number of Fisher Scoring iterations: 6Null deviance: 2.0308 on 99 degrees of freedom Residual deviance: 1.3156 on 97 degrees of freedom AIC: 798.4
deviance_inv_blocks = 1-pchisq(2.0308-1.3156, 99-97); deviance_inv_blocks [1] 0.6994No Rechazo de la Hipótesis Nula, dado que el p-valor es relativamente alto (0.6994), esto sugiere que no hay evidencia estadística significativa para rechazar la hipótesis
Prueba de Coeficientes Estimados
summary(modelo_inv_blocks)
Call:
glm(formula = Time ~ ., family = inverse.gaussian(link = "log"),
data = blocks)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.5408 0.2380 14.87 < 2e-16 ***
Child 0.1198 0.1539 0.78 0.438
Number -0.8958 0.1913 -4.68 9.5e-06 ***
Trial -0.0389 0.1286 -0.30 0.763
Shape 0.1765 0.1497 1.18 0.241
Age 0.4055 0.1557 2.60 0.011 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for inverse.gaussian family taken to be 0.01702)
Null deviance: 2.0308 on 99 degrees of freedom
Residual deviance: 1.2840 on 94 degrees of freedom
AIC: 801.9
Number of Fisher Scoring iterations: 12Intercept. Estimación: 3.5408 y p-valor: < 2e-16. Este coeficiente es estadísticamente significativo (***), lo que indica que el intercepto es diferente de cero con alta certeza.
Child.Estimación: 0.1198 y p-valor: 0.438. No es significativo (sin asterisco), lo que sugiere que la variable “Child” no tiene un efecto significativo en el tiempo en este modelo.
Number. Estimación: -0.8958 y p-valor: 9.5e-06.Este coeficiente es altamente significativo (***), lo que indica que hay una relación negativa significativa entre “Number” y “Time”. Un aumento en “Number” se asocia con una disminución en “Time”.
Trial. Estimación: -0.0389 y p-valor: 0.763.No es significativo (sin asterisco), sugiriendo que “Trial” no tiene un efecto notable en el modelo.
Shape. Estimación: 0.1765 y p-valor: 0.241. No es significativo (sin asterisco), lo que indica que “Shape” no tiene un efecto estadísticamente significativo en “Time”.
Age. Estimación: 0.4055 y p-valor: 0.011. Este coeficiente es significativo (*), lo que sugiere que “Age” tiene un efecto positivo en “Time”. Es decir, a medida que “Age” aumenta, “Time” también tiende a aumentar.
CONCLUSIONES:
Las Variables significativas “Number” y “Age” son las variables que tienen un efecto significativo en el modelo. Las Variables no significativas “Child”, “Trial” y “Shape” no son significativas, lo que sugiere que podrían no ser útiles para predecir “Time” en este contexto.
Interpretacion del coeficiente del modelo
coef(modelo_inv_blocks)(Intercept) Child Number Trial Shape Age
3.54075 0.11979 -0.89576 -0.03895 0.17651 0.40552 CONCLUSIONES
Edad muestra un efecto positivo en el tiempo que les toma a los niños apilar bloques. Esto sugiere que, a medida que aumenta la edad del niño, también aumenta el tiempo que tardan en completar la tarea. Esto podría implicar que los niños mayores pueden estar tomando más tiempo para ser más cuidadosos o estratégicos al apilar.
Número tiene un coeficiente negativo y significativo, lo que indica que a medida que aumenta la cantidad de bloques que el niño puede apilar, el tiempo requerido para completar la tarea disminuye. Esto podría reflejar una mayor habilidad o eficiencia en niños que pueden manejar más bloques, quizás debido a la práctica o experiencia.
Niño, Ensayo y Forma no mostraron un efecto significativo en el tiempo. Esto sugiere que, en este modelo, el identificador del niño y las características de los ensayos y las formas de los bloques no están directamente relacionados con la variación en el tiempo que tardan los niños en apilar.
Child
exp(0.11979)[1] 1.127Un aumento de la variable Child se asocia con un aumento del 12.7% en la variable Time
Shape
exp(0.17651)[1] 1.193Al cambiar la forma de los bloques, se espera un aumento del 19.3% en la variable Time
Age
exp(0.40552)[1] 1.5Por cada año adicional en la edad del niño, se espera que el tiempo en segundos aumente en un 50%.
Number
1/exp(-0.89576)[1] 2.449Para cada unidad adicional en el número de bloques que un niño puede apilar, el tiempo disminuye a aproximadamente 2.449 veces, lo que implica que un mayor número de bloques se asocia con una reducción en el tiempo requerido.
Trial
1/exp(-0.03895)[1] 1.04Los efectos del ensayo en el tiempo son relativamente pequeños, ya que un cambio en esta variable se asocia con un aumento del 4% en el tiempo, lo que sugiere que el ensayo no tiene un impacto considerable en el resultado.
summary(mejor_modelo_inv_blocks)
Call:
glm(formula = Time ~ Number + Age, family = inverse.gaussian(link = "log"),
data = blocks)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.746 0.163 22.96 < 2e-16 ***
Number -0.970 0.175 -5.54 2.6e-07 ***
Age 0.342 0.134 2.55 0.012 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for inverse.gaussian family taken to be 0.01808)
Null deviance: 2.0308 on 99 degrees of freedom
Residual deviance: 1.3156 on 97 degrees of freedom
AIC: 798.4
Number of Fisher Scoring iterations: 6coef(mejor_modelo_inv_blocks)(Intercept) Number Age
3.7460 -0.9703 0.3422 Number
1/exp(-0.9703)[1] 2.639Age
exp(0.3422)[1] 1.408Intercepto. Este valor indica que cuando todas las variables independientes son cero, el tiempo estimado para apilar bloques es de aproximadamente 3.746 segundos.
Número (Number). Este coeficiente es negativo, lo que sugiere que por cada bloque adicional que el niño puede apilar, el tiempo para completar la tarea disminuye. Específicamente, por cada unidad adicional en “Number”, se espera que el tiempo se reduzca a aproximadamente 2.639 veces en relación a la base. Esto resalta la eficiencia de los niños al apilar más bloques.
Edad (Age). Este coeficiente es positivo, indicando que por cada año adicional en la edad del niño, el tiempo para apilar bloques aumenta en un 40.8% (1.408 veces) manteniendo constantes las otras variables. Esto puede implicar que a medida que los niños crecen, tienden a ser más cuidadosos o a seguir un proceso más detallado al apilar.
Resumen:
El análisis indica que “Number” tiene un efecto inverso significativo sobre el tiempo, mientras que “Age” tiene un efecto directo positivo. Esto sugiere que la habilidad de apilar bloques mejora con la práctica (más bloques apilados) y que la edad también influye en el rendimiento, posiblemente debido a la madurez y la atención.
Estos hallazgos son útiles para entender cómo las habilidades de apilamiento en niños pueden desarrollarse. Podrían informar la planificación de actividades educativas que fomenten estas habilidades y considerar el impacto de la edad en el diseño de tareas similares.
(1b) Realice una estimación con: Number = 10, Trial = 2 Shape = Cube y Age = 5.
# Crear un nuevo dataframe con las variables que quieres predecir
nuevo_dato <- data.frame(
Number = 1, # Recuerde que el modelo está en términos de la variable dummy, así que si Number > 7, debe ser 1
Trial = 1, # En el modelo se asignó 1 para Trial = 2
Shape = 1, # Se asignó 1 para "Cube"
Age = 1 # En el modelo, Age también se asigna como 1, ya que está dentro del rango definido
)
# Realizar la predicción utilizando el modelo ajustado
prediccion <- predict(mejor_modelo_inv_blocks, newdata = nuevo_dato, type = "response"); prediccion 1
22.6 Según el modelo de la regresión inversamente gaussiana que se ha realizado, se estima que un niño de 5 años tardará aproximadamente 22.6 segundos en apilar 10 bloques cúbicos en la segunda prueba.
(1c) Existe datos atípicos o influyentes?. Muestre los residuales y las medidas de influencia.
Verificar datos atípicos
par(mfrow=c(1,2))
plot(abs(residuals(mejor_modelo_inv_blocks)))
abline(h=2,col = 'red')
plot(abs(residuals(mejor_modelo_inv_blocks,type = 'pearson')))
abline(h=2,col = 'red')
# Crear un dataframe con los residuos
residuos_inv_blocks = data.frame(
deviance = abs(residuals(mejor_modelo_inv_blocks)),
pearson = abs(residuals(mejor_modelo_inv_blocks, type='pearson'))
)
# Filtrar los residuos que cumplen la condición
resultados_filt_inv_blocks = residuos_inv_blocks[residuos_inv_blocks$deviance > 2 & residuos_inv_blocks$pearson > 2, ]; resultados_filt_inv_blocks[1] deviance pearson
<0 rows> (o 0- extensión row.names)La filtración de residuos en resultados_filt_inv_blocks no devolvió filas, lo que indica que no hay residuos que cumplan con las condiciones especificadas (deviance > 2 y pearson > 2). Esto sugiere que, en tu modelo, los residuos están dentro de un rango que no se considera problemático según esos umbrales.
Datos Influyentes
par(mfrow = c(1,1))
library(car)
influencePlot(mejor_modelo_inv_blocks)
StudRes Hat CookD
1 -0.4612 0.03476 0.001382
2 -1.0962 0.03476 0.004938
6 2.2521 0.03476 0.166613
24 2.9306 0.02422 0.184679
100 -2.6649 0.02422 0.011198Aunque hay algunos puntos con residuos estudentizados relativamente altos (2.2521 y 2.9306), el bajo apalancamiento y las bajas distancias de Cook sugieren que estos puntos no se pueden considerar influyentes a la luz de todas las medidas de influencia.