TAREAMODELOLINEALES OCTAVIO
OCTAVIO HERNANDEZ SATURNINO
2023-11-05
# Cargamos las librerías
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.2.3## Warning: package 'ggplot2' was built under R version 4.2.3## Warning: package 'tibble' was built under R version 4.2.3## Warning: package 'dplyr' was built under R version 4.2.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(forcats)
library(broom)
library(reshape2)##
## Attaching package: 'reshape2'
##
## The following object is masked from 'package:tidyr':
##
## smithslibrary(lmtest)## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(mgcv)## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
##
## This is mgcv 1.8-41. For overview type 'help("mgcv-package")'.library(MASS)##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## selectlibrary(modelr)##
## Attaching package: 'modelr'
##
## The following object is masked from 'package:broom':
##
## bootstrap-Caso 1 (RLS)
Se trata de determinar la pérdida de calor sufrida por cierto compuesto cuando es sometido a altas temperaturas. Los datos recogidos son los siguientes:
temperatura <- c(460, 450, 440, 430, 420, 410, 450, 440, 430,
420, 410, 400, 420, 410, 400)
perdida <- c(0.3, 0.3, 0.4, 0.4, 0.6, 0.5, 0.5, 0.6, 0.6, 0.6,
0.7, 0.6, 0.6, 0.6, 0.6)
ejer02 <- data.frame(temperatura,perdida)
ejer02En este caso ambas variables son de tipo numérico y debemos plantear un modelo de regresión lineal simple. Comenzamos con la representación gráfica de los datos.
¿Qué tipo de tendencia se aprecia en el gráfico? ¿Podemos decir que la pérdida de calor disminuye cuando aumenta la temperatura? ¿La tendencia lineal parece apropiada para este tipo de datos?
ggplot(ejer02,aes(temperatura,perdida)) +
geom_point() +
labs(x = "Temperatura", y = "Pérdida de calor") +
theme_bw() +
geom_smooth(method = "lm", se = FALSE)## `geom_smooth()` using formula = 'y ~ x'
-CONCLUSION en la grafica se aprecia que no hay normalidad Entre las
variables, puesto que los puntos estan alejados a la linea, entonces no
hay correlacion entre temperatura y perdida de calor. que perdida de
calor no depende de la temperatura.
-Estimación del modelo
Dado que sólo tenemos una predictora, la selección del modelo consiste únicamente en la valoración de la bondad del ajuste proporcionado por la variable predictora.
¿El modelo obtenido es adecuado para explicar el comportamiento de la pérdida de calor a través de la temperatura?No, ya que no hay relacion entre variables,la temperartura no dependen el comportamiento de la perdida de calor. ¿Cuál es la ecuación del modelo ajustado? calor= 2.536-0.004(Temperatura) ¿Cómo interpretamos la interceptación del modelo?
modelo de regresión lineal es:
perdida = 2.536-0.004Temperatura si la perdida incrementa,disminuye, una pm, entonces la temperatura dismuye ,incrementa. ¿Y la pendiente? Si aumentamos 100 grados ¿qué ocurre con la pérdida de calor?
perdida = 2.536-0.004(100)Temperatura ,se disminuye la perdida de calor.
#########################################
# Ajuste del modelo
ajuste <- lm(perdida ~ temperatura,data = ejer02)
# Inferencia sobre los coeficientes del modelo
tidy(ajuste)# Bondad del ajuste
glance(ajuste)# Alternativamente podríamos utilizar summary(ajuste)
########################33########3-Diagnóstico del modelo Realizamos el diagnóstico del modelo para verificar si se cumplen las hipótesis del modelo de regresión lineal simple.
El modelo ajustadono se cumple los supuestos de normalidad , la grafica se ve claro que los punto se aleja a la linea.se concluye que no hay normalidad.
Respecto del análisis gráfico ¿qué podemos decir de los gráficos obtenidos? El modelo ajustadono se cumple los supuestos de normalidad , la grafica se ve claro que los punto se aleja a la linea.se concluye que no hay normalidad.
# Valores de diagnóstico
fitinf <- fortify(ajuste)
# Análisis gráfico
diagnosis <- fitinf[,c(".fitted","temperatura",".stdresid")]
datacomp = melt(diagnosis, id.vars='.stdresid')
# Gráfico residuos, ajustados y regresora
ggplot(datacomp) +
geom_jitter(aes(value,.stdresid, colour=variable),) +
geom_smooth(aes(value,.stdresid, colour=variable),
method= mgcv::gam, se=FALSE) +
facet_wrap(~variable, scales="free_x") +
labs(x = "", y = "Residuos") +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(sample=.stdresid)) +
stat_qq() +
geom_abline() +
theme_bw()
Respecto de los tests diagnósticos ¿qué podemos decir sobre las hipótesis del modelo? El modelo ajustado para perdida de calor de se cumple los supuestos de normalidad d y suspuestos de homogeneidad de varianzas studentized Breusch-Pagan test=p-value = 0.1299,Shapiro-Wilk normality test=p-value = 0.207, independencia = 0.05 porque todos los valores p son mayores a 0.05.
¿El modelo obtened es adecuado para predecir el comportamiento de la pérdida de calor? si , por que si se cumple los supuesto de normalidad.
# Varianza constante
bptest(ajuste)##
## studentized Breusch-Pagan test
##
## data: ajuste
## BP = 2.294, df = 1, p-value = 0.1299# Normalidad
shapiro.test(fitinf$.stdresid)##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid
## W = 0.92206, p-value = 0.207# Independencia
dwtest(ajuste)##
## Durbin-Watson test
##
## data: ajuste
## DW = 1.3434, p-value = 0.05764
## alternative hypothesis: true autocorrelation is greater than 0-Predicción del modelo Pasamos a la predicción de la media para el modelo ajustado.
¿Cuál es estimación de la pérdida de calor para las temperaturas de 400, 420, 440, y 460? de 400 a 460 se dismuyen.
¿Cuál es el rango de confianza para la predicción en cada una de las temperaturas anteriores? 400,460 ¿Resulta posible calcular la pérdida de calor para una temperatura de 480? la grafica tiene forma de dismucion , entonces en 480 esta en forma de disminución ¿Cuál sería el intervalo de confianza para dicha predicción? 400,480
# Predicción
newdata <- data.frame(temperatura = seq(min(ejer02$temperatura),
max(ejer02$temperatura), .01))
# Obtenemos la predicción para el modelo ajustado
newdata <- data.frame(newdata, predict(ajuste, newdata,
interval="confidence"))
# Gráfico de la predicción
ggplot(newdata, aes(x = temperatura, y = fit)) +
geom_line() +
geom_ribbon(aes(ymax = upr, ymin = lwr), alpha = 1/5) +
labs(x = "Temperatura",
y = "Pérdida de calor",
title = "Bandas de predicción") +
theme_bw()
-\(Caso 2 (RLM)\) Los datos muestran el porcentaje de calorías totales obtenidas de carbohidratos complejos, para veinte diabéticos dependientes de insulina que habían seguido una dieta alta en carbohidratos durante seis meses. Se consideró que el cumplimiento del régimen estaba relacionado con la edad (en años), age, el peso corporal (relativo al peso “ideal” para la altura), weight, y otros componentes de la dieta como el porcentaje de proteínas ingeridas. Los datos corresponden con la tabla 6.3 de Dobson (2002). u
# Lectura de datos
ejer03 <- read_csv("https://goo.gl/Grm8xM", col_types = "dddd")En este caso todas las variables son de tipo numérico, y debemos plantear un modelo de regresión lineal múltiple donde tratamos de explicar el comportamiento de carbohydrate en función del conjunto de posibles predictoras. Comenzamos con la representación gráfica de los datos a través de gráficos de dispersión individualizados de la respuesta con cada posible predictora.
¿Qué variable o variables parecen influir más, de forma individual, en el porcentaje de calorías totales obtenidas de carbohidratos complejos? PROTEINAS ¿Qué tipo de tendencia se observa para cada una de las predictoras? la grafica en años va disminuyendo lento. la grafica en weingt va disminuyendo de forma rapida. la grafica en proteinas va incrementado rapido. ¿Son tendencia lineales o podríamos sospechar la presencia de tendencias no lineales? si son tendencias lineales.porque tiene comportamiento de tendencia de normalidad.
# Gráfico inicial
datacomp = melt(ejer03, id.vars='carbohydrate')
# Gráfico respecto de cada predictora
ggplot(datacomp) +
geom_jitter(aes(value,carbohydrate, colour=variable),) +
geom_smooth(aes(value,carbohydrate, colour=variable),
method= mgcv::gam, se=FALSE) +
facet_wrap(~variable, scales="free_x") +
labs(x = "", y = "carbohydrate") +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
-Estimación del modelo
En primer lugar procedemos con la selección del modelo. Posteriormente haremos el análisis inferencial de dicho modelo. Partimos del modelo completo con todas las predictoras y utilizamos el procedimiento por pasos para la selección de efectos en el modelo final.
¿Qué variables predictoras aparecen el modelo final? carbohydrate - weight + protein
Justifica la posible eliminación de alguna o algunas de las variables predictoras?
el age es una variable que no tiene correlación con carbohydrate, pues la variable age no aporta nada, y es posible eliminarlo sin afectar ala variable respuesta carbohydrate, al final el modelo se queda con las siguientes variables carbohydrate - weight + protein
# Ajuste del modelo
ajuste <- lm(carbohydrate ~ age + weight + protein,
data = ejer03)
# Selección del modelo
ajuste <- step(ajuste)## Start: AIC=74.92
## carbohydrate ~ age + weight + protein
##
## Df Sum of Sq RSS AIC
## - age 1 38.36 606.02 74.224
## <none> 567.66 74.916
## - weight 1 265.91 833.57 80.600
## - protein 1 337.34 905.00 82.244
##
## Step: AIC=74.22
## carbohydrate ~ weight + protein
##
## Df Sum of Sq RSS AIC
## <none> 606.02 74.224
## - weight 1 252.63 858.65 79.193
## - protein 1 305.40 911.42 80.385¿Consideras que el modelo obtenido es adecuado para explicar el porcentaje de calorías totales obtenidas de carbohidratos complejos? el 44% de la variabilidad de carbohidrato es explicado por el modelo de regresion lineal. Justifica tu respuesta. ¿Los efectos presentes en el modelo se pueden considerar como relevantes?
¿Cuál es el modelo ajustado? Interpreta los coeficientes de dicho modelo. Con un nivel de significancia de 𝛼 = 0.05, existe evidencia que las medias del del porcentaje carbohidratos con weingth , protein no son diferente (𝐹𝑐(1, 17) = 35.64, 𝑝 > 0.009)
# Ajuste del modelo final
ajuste <- lm(carbohydrate ~ weight + protein,
data = ejer03)
# Inferencia
tidy(ajuste)# Bondad del ajuste
glance(ajuste)# Tabla ANOVA
anova(ajuste)-Diagnóstico del modelo En este caso nos centramos únicamente en los procedimientos numéricos y tests. ¿Se detecta alguna observación influyente? El modelo ajustado para carbohidrato se cumple los supuestos de normalidad y supuestos de varianzas de homogeneidad
¿Que indican los tests de diagnóstico? hay homogeneidad de varianza,se rechaza Ho, puesto que valores p son mayores que 0.05,que Las medias de carbohidrato en weingth y proteins NO son iguales. ¿Podemos utilizar el modelo ajustado para la predicción del porcentaje de de calorías totales obtenidas de carbohidratos complejos?no, tomando en cuenta de lo anterior , entonces El modelo ajustado para porcentaje de carbohidrato se cumple los supuestos de normalidad y supuestos de varianzas de homogeneidad, independencia de autocorrelacion.
# Diagnóstico del modelo
fitinf <- fortify(ajuste)
fitinf# Análisis gráfico
diagnosis <- fitinf[,c(".fitted","protein", "weight", ".stdresid")]
datacomp = melt(diagnosis, id.vars='.stdresid')
# Tests estadísticos
# Varianza constante
bptest(ajuste)##
## studentized Breusch-Pagan test
##
## data: ajuste
## BP = 0.19392, df = 2, p-value = 0.9076# Normalidad
shapiro.test(fitinf$.stdresid)##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid
## W = 0.95111, p-value = 0.3843# Independencia
dwtest(ajuste)##
## Durbin-Watson test
##
## data: ajuste
## DW = 2.1074, p-value = 0.5579
## alternative hypothesis: true autocorrelation is greater than 0-Predicción Construimos un rango de predicción para las proteínas y trabajamos con un peso medio para obtener la gráfica de predicción.
¿Cómo interpretamos este gráfico de predicción? que ahi correlacion entre carbohidrato y proteinas, la grafica es creciente, va aumentando conforme va avanzado.
¿Tendría sentido reconvertir la variable peso en un factor ordinal y ajustar el modelo correspondiente?no tiene caso, la grafica se ve perfectamente que la variable proteinas es explicada por la variable respuesta de carbohidrato. Justifica tu respuesta. ¿Qué tipo de modelo tendríamos entonces? modelo de regresion linea en forma positivo, incrementado. ¿Sería más o menos informativo que el modelo con todas las predictoras numéricas? ¿Qué diferencias en el estudio del modelo deberíamos tener en cuenta? Ajusta un nuevo modelo considerando tres grupos en la variable peso identificados como menos de 100, entre 100 y 120, y más de 120.
newdata <-expand.grid(protein = round(seq(min(ejer03$protein),
max(ejer03$protein),
length = 10),0),
weight = round(mean(ejer03$weight),0))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(ajuste, newdata,
interval="confidence"))
# Gráfico para cada valor de weight
ggplot(newdata, aes(x = protein, y = fit)) +
geom_line() +
geom_ribbon(aes(ymax = upr, ymin = lwr), alpha = 1/5) +
labs(x = "Protein", y ="Carbohydrate",
title = "Mean Weight") +
facet_wrap(~ weight, scales="free_x") +
theme_bw()
Construcción de la nueva variable
# Creación de la nueva varaible
ejer03$weight.fac <- cut(ejer03$weight,breaks=c(50,100,120,150))-\(Caso 3 (RLS)\)
Los datos recogen la respuesta de un pasto de gramíneas y leguminosas a varias cantidades de fertilizante de fósforo (datos de D. F. Sinclair). En concreto se mide el rendimiento total (en kilogramos por hectárea), de pasto y leguminosa juntos, y la cantidad de potasio ([Math Processing Error] ) utilizada (en kilogramos por hectárea). Los datos corresponden con la tabla 6.16 de Dobson (2002).
# Lectura de datos
ejer04 <- read_csv("https://goo.gl/AOikQU", col_types = "dd")En este caso tenemos las dos variables de tipo numérico y planteamos un modelo de regresión lineal simple.
¿Cómo afecta la cantidad de potasio en el rendimiento total del pasto? en la grafica no hay normalidad, quiere decir que no hay correlacion entre las las variables yield y K ¿Podemos considerar la relación entre ambas variables de tipo lineal? que no hay relación entre las variables que yield no depende de la variable K, por lo tanto no hay relación entre variables.
# Gráfico inicial
ggplot(ejer04,aes(K,yield)) +
geom_point() +
labs(x = "K", y = "yield") +
theme_bw() +
geom_smooth(method = "lm", se = FALSE)## `geom_smooth()` using formula = 'y ~ x'
-Estimación del modelo.
Construimos y evaluamos la calidad del ajuste obtenido.
¿Podemos decir que la cantidad de potasio explica el rendimiento total del pasto? Justifica tu respuesta. El 77% de la variabilidad del porcentaje yield en k es explicado por el modelo de regresion lineal.r. ¿Que relación se ha establecido entre ambas variables? que hay buena relacion ¿Cuál sería el rendimiento del pasto si no utilizamos el potasio? yield=2644+48.36*(potacio), entonces si no utilizamo el potacio en tonces obtenemos yield =2644+48.36
¿Cómo aumenta el rendimiento por cada aumentar en 10 kilogramos por hectárea la cantidad de potasio?yield=2644+48.36*(10kiligramo)
# Ajuste del modelo
ajuste <- lm(yield ~ K,data = ejer04)
# Análsis inferencial sobre los coeficientes del modelo
tidy(ajuste)# Bondad del ajuste
glance(ajuste)-Diagnóstico del modelo Realizamos el diagnóstico del modelo tanto mediante procedimientos gráficos como los tests de diagnóstico.
¿Cómo interpretamos los gráficos de diagnóstico? El modelo ajustado para yield se cumple los supuestos de normalidad y supuestos de varianzas de homogeneidad
¿Encontramos alguna dificultad con las hipótesis del modelo?si, en la gráfica los punto se aleja de la linea, se ve que no hay normalidad.pero con el test, uno se da cuenta de que hay normalidad,
¿Podemos considerar que le modelo obtenido es adecuado para explicar la relación entre la cantidad de potasio y el rendimiento del pasto? que se cumple la normalidad la homogeneidad de varianza y independencia.
# Diagnóstico del modelo
fitinf <- fortify(ajuste)
# Análisis gráfico
diagnosis <- fitinf[,c(".fitted","K",".stdresid")]
datacomp = melt(diagnosis, id.vars='.stdresid')
# Gráfico residuos, ajustados y regresora
ggplot(datacomp) +
geom_jitter(aes(value,.stdresid, colour=variable),) +
geom_smooth(aes(value,.stdresid, colour=variable),
method= mgcv::gam, se=FALSE) +
facet_wrap(~variable, scales="free_x") +
labs(x = "", y = "Residuos") +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(sample=.stdresid)) +
stat_qq() +
geom_abline() +
theme_bw()
# Varianza constante
bptest(ajuste)##
## studentized Breusch-Pagan test
##
## data: ajuste
## BP = 0.33734, df = 1, p-value = 0.5614# Normalidad
shapiro.test(fitinf$.stdresid)##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid
## W = 0.97288, p-value = 0.6792# Independencia
dwtest(ajuste)##
## Durbin-Watson test
##
## data: ajuste
## DW = 2.1839, p-value = 0.6932
## alternative hypothesis: true autocorrelation is greater than 0-Predicción
Calculamos un grid de valores de potasio y obtenemos el correspondiente modelo de predicción. ¿Cuáles son la predicciones del rendimiento del pasto para los valores de potasio de 20, 30, y 40? de intervalo(0,50) prediciones va incrementado, es modelo lineal
newdata <- data.frame(K = seq(min(ejer04$K), max(ejer04$K), .5))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(ajuste, newdata,
interval="confidence"))
# Gráfico de la predicción
ggplot(newdata, aes(x = K, y = fit)) +
geom_line() +
geom_ribbon(aes(ymax = upr, ymin = lwr), alpha = 1/5) +
labs(x = "K", y = "Yield", title = "Bandas de predicción") +
theme_bw()
-\(Caso 4 (RLM)\) Es bien sabido que la concentración de colesterol en el suero sanguíneo aumenta con la edad, pero es menos claro si el nivel de colesterol también está asociado con el peso corporal. Los datos muestran para una treinta de mujeres el colesterol sérico (milimoles por litro), la edad (años) y el índice de masa corporal (peso dividido por la altura al cuadrado, donde el peso se midió en kilogramos y la altura en metros). Los datos corresponden con la tabla 6.17 de Dobson (2002).
Lectura de datos
# Lectura de datos
ejer05 <- read_csv("https://goo.gl/EKXWRc",
col_types = "ddd")Debemos plantear un modelo de regresión lineal múltiple donde tratamos de explicar la concentración de colesterol en el suero sanguíneo en función de la edad y del peso corporal. Realizamos gráficos de dispersión individual para cada posible predictora.
¿Qué conclusiones podemos extraer de este gráfico? grafica (age) no correlacion porque los puntos estan desordenados alejados en las lineas. grafica(bmi)no correlacion porque los puntos estan desordenados alejados en las lineas ¿Las dos predictoras pueden tener efecto sobre la concentración de colesterol? no.se tiene que hacer prueva de normalidad y independencia para saber el comportamiento de la grafica.
datacomp = melt(ejer05, id.vars='chol')
# Gráfico respecto de cada predictora
ggplot(datacomp) +
geom_jitter(aes(value,chol, colour=variable),) +
geom_smooth(aes(value,chol, colour=variable),
method= mgcv::gam, se=FALSE) +
facet_wrap(~variable, scales="free_x") +
labs(x = "", y = "Chol") +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
-Estimación del modelo Procedemos con la estimación del mejor modelo realizando un proceso de selección de variables.
¿Qué predictoras contiene el modelo resultante del proceso de selección de efectos? chol ~ age + bmi el modelo es adecuado para la expliacacion de la variable respuesta.
¿El modelo obtenido es adecuado para explicar la concentración de colesterol? Justifica tu respuesta. si.El 46% de la variabilidad es explicado pòr el chol con age+ bmil es explicado pr el modelo de regresion lineal múltiple ¿Cuál es la ecuación de modelo ajustado?
chol=-0.7398262 + 0.04096772*(age)+0.20137026(bmi) ¿Cómo aumenta la concentración del colesterol para dos mujeres con un diferencia de edad de cinco años?
¿y de un volumen de masa corporal de de cinco unidades? ¿Tiene sentido una interceptación negativa para este modelo?
¿Cuál sería la posible explicación? Vamos a probar a ajustar un modelo sin interceptación y compararemos los resultados con el modelo con interceptación.
llegamos ala conclusion de que no se rechaza ho.Con un nivel de significancia de 𝛼 = 0.05 Hay efecto de age y bmi en chol.por p_valor menor que 0.05.
# Ajuste del modelo
ajuste <- lm(chol ~ age + bmi ,
data = ejer05)
# Selección del modelo
ajuste <- step(ajuste)## Start: AIC=2.36
## chol ~ age + bmi
##
## Df Sum of Sq RSS AIC
## <none> 26.571 2.3583
## - bmi 1 5.0655 31.636 5.5931
## - age 1 8.8907 35.461 9.0173# Inferencia
tidy(ajuste)# Bondad del ajuste
glance(ajuste)Modelo sin interceptación.
¿Resulta más adecuado este modelo? Justifica tu respuesta y responde a las preguntas sobre el modelo planteadas anteriormente. a lo que llagamos de acuerdo el modelo presentado(chol ~ age + bmi - 1) llegamos a la conclusion que el 97% de la variabilidad es explicada por el modelo 1(chol ~ age + bmi - 1) chol= 0.04135672(age)+ 0.16889397(bmi)-1 Con un nivel de significancia de 𝛼 = 0.05 Hay efecto de age y bmi en chol.por p_valor menor que 0.05,y no se rechaza ho.
# Ajuste del modelo sin interceptación
ajuste2 <- lm(chol ~ age + bmi - 1,
data = ejer05)
# Selección del modelo
ajuste2 <- step(ajuste2)## Start: AIC=0.53
## chol ~ age + bmi - 1
##
## Df Sum of Sq RSS AIC
## <none> 26.720 0.5269
## - age 1 9.109 35.829 7.3272
## - bmi 1 29.619 56.340 20.9062# Inferencia
tidy(ajuste2)# Bondad del ajuste
glance(ajuste2)# Comparación de ambos modelos
anova(ajuste,ajuste2)ajuste <- ajuste2-Diagnóstico del modelo Realizamos el diagnóstico del modelo tanto mediante procedimientos gráficos como los tests de diagnóstico.
¿Que opinas sobre los gráficos de diagnóstico? que los tres gráficos tiene un comportamiento similar. ¿Podríamos considerar que existen observaciones influyentes? Realizamos el análisis correspondiente. de acuerdo de la gráfica podemos ver que hay normalidad.que los residuos tiene distribución normal. ¿Que conclusión extraemos de este análisis? ¿Debemos considerar la eliminación de alguna observación? Justifica tu respuesta.
# Diagnóstico del modelo
fitinf <- fortify(ajuste)
# Análisis gráfico
diagnosis <- fitinf[,c(".fitted","age", "bmi", ".stdresid")]
datacomp = melt(diagnosis, id.vars='.stdresid')
# Gráfico residuos, ajustados y regresora
ggplot(datacomp) +
geom_jitter(aes(value,.stdresid, colour=variable),) +
geom_smooth(aes(value,.stdresid, colour=variable),
method= mgcv::gam, se=FALSE) +
facet_wrap(~variable, scales="free_x") +
labs(x = "", y = "Residuos") +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(sample=.stdresid)) +
stat_qq() +
geom_abline() +
theme_bw()
# Análisis de influencia
fitinf$.cooksd > 1## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSE-Tests diagnósticos.
¿Podemos considerar como adecuado el modelo obtenido? Justifica tu respuesta. se cumple los supuestos de normalidad(shapiro_wilk =0.20) y supuestos de varianzas de homogeneidad(studentized Breusch_Pagan test=0.059)y por ultimo IDEPENDENCIA ( Durbin_Watson= 0.89). Hay correlaciones entre las variables.
# Varianza constante
bptest(ajuste)##
## studentized Breusch-Pagan test
##
## data: ajuste
## BP = 3.6464, df = 1, p-value = 0.05619# Normalidad
shapiro.test(fitinf$.stdresid)##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid
## W = 0.95325, p-value = 0.2064# Independencia
dwtest(ajuste)##
## Durbin-Watson test
##
## data: ajuste
## DW = 2.4589, p-value = 0.8975
## alternative hypothesis: true autocorrelation is greater than 0-Predicción Construimos un grid de valores para la edad y trabajamos con el valor mínimo, máximo y la media del índice de masa corporal para representar la predicción asociada.
¿Cómo podemos interpretar los modelos de predicción obtenidos? La grafica tiene el mismo comportamiento ¿Cuál es la diferencia de predicción para una edad de 50 años? la grafica de 19 son 5 la grafica de 23 son 6 lagrafica de 29 son 7 la diferencia de los tres graficas es uno.
# Predicción
newdata <-expand.grid(age = round(seq(min(ejer05$age),
max(ejer05$age),
length = 10),0),
bmi = round(c(min(ejer05$bmi),
mean(ejer05$bmi),
max(ejer05$bmi)),0))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(ajuste, newdata,
interval="confidence"))
# Gráfico para cada valor
ggplot(newdata, aes(x = age, y = fit)) +
geom_line() +
geom_ribbon(aes(ymax = upr, ymin = lwr), alpha = 1/5) +
labs(x = "Age", y ="Cholesterol", title = "BMI") +
facet_wrap(~ bmi, scales="free_x") +
theme_bw()
-\(Caso 5 (RLM)\) En un estudio sobre las diferentes clases de queso cheddar que se fabrican en LaTrobe Valley de Victoria, Australia, se analizaron muestras de queso por su composición química: concentración de ácido acético (escala logarítmica); concentración de sulfuro de hidrógeno (escala logarítmica); y la concentración de ácido láctico. Por otro lado, se paso una muestra de cada uno de ellos a un conjunto de catadores y se registro la puntuación obtenida por cada uno de ellos. Estamos interesados en relacionar la puntuación final de los catadores con los resultados del análisis químico.
# Lectura de datos
ejer06 <- read_csv("https://goo.gl/V4lDNs", col_types = "dddd")Nos encontramos ante un modelo de regresión lineal múltiple (todas las variables son de tipo numérico) donde tratamos de explicar la puntuación final de cada queso en función de la concentración de ácido acético, concentración de sulfuro de hidrógeno, y la concentración de ácido láctico. cargamos los datos y realizamos los gráficos de dispersión individuales.
¿Cómo podemos interpretar los gráficos obtenidos? los tres tienen comportamientos similares.y son modelos lineales generalizados.por ser asi.
¿tenemos alguna indicación de las variables que pueden resultar más relevantes para explicar la puntuación final obtenida?los tres son modelos lineales, entonces los tres tiene relacion en comun.
# Gráfico inicial
datacomp = melt(ejer06, id.vars='Taste')
# Gráfico respecto de cad predictora
ggplot(datacomp) +
geom_jitter(aes(value,Taste, colour=variable),) +
geom_smooth(aes(value,Taste, colour=variable),
method= mgcv::gam, se=FALSE) +
facet_wrap(~variable, scales="free_x") +
labs(x = "", y = "Taste") +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
-Estimación del modelo
Construimos el mejor modelo posible en función de las variables
consideradas. ¿Cuáles son las predictoras involucradas en el modelo
final? Traste=H2S + Lactic H2S Y Lactic son las variables que aportan
mas a trate la caetic es la variable que no aporta nada, y al final solo
se reduce, se quita la variable que no aporta nada, y se obtine
Traste=H2S + Lactic ¿Cuál es la capacidad explicativa de dicho modelo?
Traste=-27.59+3.94(H2S)+19.88(Lactic) Si elH2s incrementa
(disminuye) una ppm, entonces Traste disminuye (incrementa)en promedio
3.9 Si Lactic incrementa (disminuye) una ppm, entonces Traste disminuye
(incrementa)en promedio 19.88 ¿Cuáles son los coeficientes del modelo
obtenido? El 65% de la variabilidad Traste es explicado en H2S Y lactic
por el modelo de regresion lineaL multiple ¿Tienen sentido o deberíamos
modificar nuestro modelo? el modelo tiene sentido.El 65% de la
variabilidad Traste es explicado en H2S Y lactic por el modelo de
regresion lineaL multiple ¿Cuál es tu propuesta de mejora? el modelo
esta mejorado, no se tiene que hacer nada, ni eliminar variables , otra
recomendacion podemos hacer la transformacion Box Cox, o con
logaritmo.
# Ajuste del modelo
ajuste <- lm(Taste ~ Acetic + H2S + Lactic,
data = ejer06)
# Selección del modelo
ajuste <- step(ajuste)## Start: AIC=142.64
## Taste ~ Acetic + H2S + Lactic
##
## Df Sum of Sq RSS AIC
## - Acetic 1 0.55 2669.0 140.65
## <none> 2668.4 142.64
## - Lactic 1 533.32 3201.7 146.11
## - H2S 1 1007.66 3676.1 150.25
##
## Step: AIC=140.65
## Taste ~ H2S + Lactic
##
## Df Sum of Sq RSS AIC
## <none> 2669.0 140.65
## - Lactic 1 617.18 3286.1 144.89
## - H2S 1 1193.52 3862.5 149.74# Inferencia
tidy(ajuste)# Bondad del ajuste
glance(ajuste)Construimos el nuevo modelo y lo comparamos con el anterior. ¿Podemos considerarlo mejor que el anterior? el modelo Taste = acetic +H2S + Lactic-1 lo interpretamos en ecuacion Taste= -5.45Acetic +4.5H2S+ 19.12*Lactic ¿Cómo interpretamos el comportamiento de cada predictora en este nuevo modelo? 88% de la variabilidad Taste es explicado por el diseño de modelo lineal multiple. Con un nivel de significancia de 𝛼 = 0.05, se rechaza Ho,porque p_valor es menor que 0.05 hay efecto de Acetic ,H2S, Lactic-1 en Taste. ,
# Ajuste del modelo
ajuste2 <- lm(Taste ~ Acetic + H2S + Lactic - 1,
data = ejer06)
# Selección del modelo
ajuste2 <- step(ajuste2)## Start: AIC=143.02
## Taste ~ Acetic + H2S + Lactic - 1
##
## Df Sum of Sq RSS AIC
## <none> 2888.1 143.01
## - Lactic 1 505.18 3393.3 145.85
## - Acetic 1 713.67 3601.8 147.64
## - H2S 1 1589.14 4477.3 154.17# Inferencia
tidy(ajuste2)# Bondad del ajuste
glance(ajuste2)-Diagnóstico del modelo Procedemos con el diagnóstico del último modelo ajustado. Comenzamos con el análisis gráfico y de influencia.
¿Qué podemos concluir del análisis realizado? por lo anterior concluimos ,que el modelo ajustado cumple con la normalidad.
ajuste <- ajuste2
# Diagnóstico del modelo
fitinf <- fortify(ajuste)
# Análisis gráfico
diagnosis <- fitinf[,c(".fitted", "H2S", "Lactic", "Acetic", ".stdresid")]
datacomp = melt(diagnosis, id.vars='.stdresid')
# Gráfico residuos, ajustados y regresora
ggplot(datacomp) +
geom_jitter(aes(value,.stdresid, colour=variable),) +
geom_smooth(aes(value,.stdresid, colour=variable),
method= mgcv::gam, se=FALSE) +
facet_wrap(~variable, scales="free_x") +
labs(x = "", y = "Residuos") +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(sample=.stdresid)) +
stat_qq() +
geom_abline() +
theme_bw()
# Análisis de influencia
fitinf$.cooksd > 1## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSERealizamos los tests diagnósticos.
¿Es válido el modelo construido? Justifica tu respuesta. El modelo ajustado para Taste se cumple los supuestos de normalidad(shapiro-wikl =0.17) se cumple homogeneidad de varianza(studentized Breusch=0.75)Independenica(Darvin-Watson=0.12)
# Varianza constante
bptest(ajuste)##
## studentized Breusch-Pagan test
##
## data: ajuste
## BP = 0.57547, df = 2, p-value = 0.75# Normalidad
shapiro.test(fitinf$.stdresid)##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid
## W = 0.95033, p-value = 0.1726# Independencia
dwtest(ajuste)##
## Durbin-Watson test
##
## data: ajuste
## DW = 1.6179, p-value = 0.1202
## alternative hypothesis: true autocorrelation is greater than 0-Predicción Construimos un grid de valores para el ácido láctico y tomamos los valores medios para las otras dos predictoras.
¿Cuál es la predicción de la puntuación para para los valores medios de las tres predictoras? Si el Acetic incrementa (disminuye) una ppm, entonces Taste disminuye (incrementa) en promedio -5.45 Si el H2S incrementa (disminuye) una ppm, entonces Taste disminuye (incrementa) en promedio 4.2 Si el Lactic incrementa (disminuye) una ppm, entonces Taste disminuye (incrementa) en promedio 19.12 ¿Qué variable crees que influirá más en la predicción final? la variable Lactic ¿Cómo podríamos comprobar este hecho? nose necesita cjhecar solo con ver predictores ,es la variable que aporta mas.
newdata <-expand.grid(H2S = mean(ejer06$H2S),
Lactic = seq(min(ejer06$Lactic),
max(ejer06$Lactic),
by = 0.2),
Acetic = mean(ejer06$Acetic))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(ajuste, newdata,
interval="confidence"))
# Gráfico para cada valor de weight
ggplot(newdata, aes(x = Lactic, y = fit)) +
geom_line() +
geom_ribbon(aes(ymax = upr, ymin = lwr), alpha = 1/5) +
labs(x = "Lactic", y ="Taste") +
theme_bw()
-\(Caso 6 (ANOVA 1via)\) En una investigación para la reducción de peso corporal se registran los pesos, en kilogramos, de veinte hombres antes y después de la participación en un programa de “pérdida de cintura”. Los datos corresponden con la tabla 2.8 de Dobson (2002). Para analizar adecuadamente estos datos juntamos las columnas correspondientes a antes y después en una columna con un factor que indica a cuando corresponde la observación.
# Lectura de datos
previo <- read_csv("https://goo.gl/jWGurk", col_types = "idd")
# Para construir modelos trasformamos los datos originales en la forma habitual
ejer07 <- previo %>% gather(`before`,`after`,
key = "Time",
value = "Weight")Representamos los datos en la forma habitual.
¿Qué conclusión podemos extraer de este gráfico? que las mediana son casi similares after y before ¿Qué resultado esperas que proporcione el modelo que puedas ajustar? que hay variabilidad entre las variables, puede explicar a la variable respùesta
# Gráfico de cajas
ggplot(ejer07,aes(x = Time, y = Weight)) +
geom_boxplot() +
theme_bw()
Estimación del modelo Atendiendo a la tabla ANOVA resultante para este modelo ¿Cuál es la conclusión del modelo obtenido? (𝐹𝑐(1, 38) = 0.41 , 𝑝 < 0.5). time no es significativo (F1,13 = 90,1, p < 0,001) en el modelo de regresion lineal simple. Eso quiere decir queel time no permite explicar weight
¿Es necesario seguir trabajando sobre este modelo? Justifica tu respuesta. no, se puede agregar otro variable o hacer una transformacion.
# Ajuste del modelo
ajuste <- lm(Weight ~ Time, data = ejer07)
# Tabla ANOVA
anova(ajuste)# Análsis inferencial sobre los coeficientes del modelo
tidy(ajuste)-\(Caso 7 (ANOVA 1via)\)
Se realiza una investigación para conocer los niveles de fosfato inorgánico en plasma (mg / dl) una hora después de una prueba de tolerancia a la glucosa estándar para sujetos obesos, con o sin hiperinsulinemia (H-O, y N-O respectivamente), y controles (C). Los datos corresponden con la tabla 6.18 de Dobson (2002).
# Lectura de datos
ejer08 <- read_csv("https://goo.gl/3L4EtK", col_types = "cd")En este caso tenemos un modelo ANOVA de un vía (grupo de pacientes) para explicar el nivel de fosfato inorgánico en plasma. Cargamos los datos y representamos de la forma habitual.
¿Qué conclusión podemos extraer de este gráfico? que los tres graficos tiene mediana diferente.H-O tiene mediana alta,y C tiene mediana baja
¿Consideras adecuado plantear un modelo para tratar de conocer si el nivel medio de fosfato es diferente entre los grupos considerados?si.
# Gráfico de cajas
ggplot(ejer08,aes(x = Group, y = phosphate)) +
geom_boxplot() +
theme_bw()
-Estimación del modelo Construimos el modelo correspondiente a esta situación.
¿Qué conclusión podemos extraer de la tabla ANOVA obtenida?
Group es significativo(𝐹𝑐(2, 28) = 11.56 , 𝑝 < 0.0002). en el modelo de regresion lineal simple. Eso quiere decir queel Group permite explicar Fosfato.
¿Cuál es la media estimada del nivel de fosfato para cada grupo? fosfato = 2.78+ GroupH-O1.16 + GroupN-O0.65 Si el GroupH-O incrementa (disminuye) una ppm, entonces Fosfato disminuye (incrementa) en promedio 1.16 Si el GroupN-O incrementa (disminuye) una ppm, entonces Fosfato disminuye (incrementa) en promedio 0.65
¿Podemos considerar que las medias en todos los grupos son diferentes? si.
# Ajuste del modelo
ajuste <- lm(phosphate ~ Group, data = ejer08)
# Tabla ANOVA
anova(ajuste)# Análsis inferencial sobre los coeficientes del modelo
tidy(ajuste)Realizamos el test de comparación múltiples.
¿Qué conclusiones podemos extraer del análisis realizado? H-O es estadisticamente igual a c mientras n-o es diferente a c ¿tendrá interés el diseño experimental llevado a cabo? Justifica tu respuesta si.
# Comparaciones múltiples
TukeyHSD(aov(phosphate ~ Group, data = ejer08))## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = phosphate ~ Group, data = ejer08)
##
## $Group
## diff lwr upr p adj
## H-O-C 1.1621212 0.5642291555 1.7600133 0.0001347
## N-O-C 0.6541667 0.0003963287 1.3079370 0.0498380
## N-O-H-O -0.5079545 -1.1735054769 0.1575964 0.1607135Diagnóstico del modelo En este caso hemos de tener en cuenta que al tratarse de un modelo ANOVA la verificación de las hipótesis son un poco más laboriosas. Nos centramos únicamente en los tests diagnósticos.
¿Qué podemos concluir del proceso de diagnóstico del modelo? su varianza de homogeneidad Bartlett test of homogeneity of variances, No existe evidencia con un α=.05 que indique que las varianzas de los phospathe en group sean diferentes.
¿Podemos utilizar el modelo construido para estimar la media del nivel de fosfato en función del grupo al que pertenecen los sujetos? el modelo ajustado se cumple el supusto de normalidad.por lo tanto en homegeneidad de varianza se rechaza ho. y se llega la conclusion. No existe evidencia con un α=.05 que indique que las varianzas de los phospathe en group sean diferentes.
# Diagnóstico
fitinf <- fortify(ajuste)
# Varianza constante
bartlett.test(phosphate ~ Group, data = ejer08)##
## Bartlett test of homogeneity of variances
##
## data: phosphate by Group
## Bartlett's K-squared = 4.663, df = 2, p-value = 0.09715# Normalidad (Obtenemos los pvalores del test de normalidad)
shapiro.test(fitinf$.stdresid[fitinf$Group == "C"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Group == "C"]
## W = 0.96, p-value = 0.7838shapiro.test(fitinf$.stdresid[fitinf$Group == "H-O"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Group == "H-O"]
## W = 0.95486, p-value = 0.7064shapiro.test(fitinf$.stdresid[fitinf$Group == "N-O"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Group == "N-O"]
## W = 0.88368, p-value = 0.2041-Predicción En esta caso la predicción se basa únicamente en el cálculo de los intervalos de confianza de la media para cada uno de los grupos considerados.
¿Qué podemos decir de los intervalos de predicción obtenidos?
¿Podríamos utilizar esta predicción como diagnóstico médico? Justifica tu respuesta.
# Predicción
# Secuencia de valores de predicción
newdata <- data.frame(Group = unique(ejer08$Group))
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata,
interval="confidence"))
# Gráfico resultante (predicción y bandas de predicción).
# Ordenamos las clases de acuerdo al valor de fit
ggplot(newdata, aes(x = fct_reorder(Group,fit), y= fit)) +
geom_errorbar(aes(ymin=lwr, ymax=upr), width=.1) +
geom_line() +
geom_point() +
labs(x = "Grupo", y = "Nivel de fósfato") +
coord_flip() +
theme_bw()## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
-\(Caso 8 (ANCOVA)\) Se lleva a cabo un estudio sobre el contenido promedio de grasa (en porcentaje) en la leche del ganado de cinco razas distintas canadienses. Para ello se consideran veinte ejemplares de pura raza (diez de dos años y diez maduras de más de cuatro años de cada una de cinco razas.
ejer09 <- read_csv("https://goo.gl/J2ZKWK", col_types = "dcc")Nos encontramos ante un diseño experimental con dos posibles variables predictoras de tipo categórico. Debemos plantear un modelo ANOVA con dos vías de clasificación. Recordemos que en estas situaciones el estudio principal se centra en primer lugar en determinar si existe interacción entre los factores considerados, es decir, si la media del contenido de grasa en la leche depende simultáneamente de la raza de la oveja y de su edad. Cargamos los datos y realizamos el gráfico de cajas conjunto y el de interacción de las medias de cada combinación de ambos factores.
¿Qué conclusiones podemos extraer del gráfico de cajas conjunto? se puede apreciar que bread es mas efectivo para cada age. ¿Qué consideras más relevante del gráfico obtenido? ¿Qué podemos decir a la vista del gráfico de interacción? que tiene el mismo comportamiento desde inicio hasta final.resulta ser medias iguales. ¿Qué tipo de modelo parecería más adecuado a la vista de estos datos?
butterfart - bread
# Diagrama de cajas
ggplot(ejer09,aes(x = Bread, y = Butterfat, color = Age)) +
geom_boxplot() +
theme_bw()
# Gráfico de interacción de medias
ggplot(ejer09, aes(x = Bread, y = Butterfat,
group = Age,
color = Age)) +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.y = mean, geom = "line") +
theme_bw()## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
-Estimación del modelo Procedemos con la estimación del modelo completo considerando la posible interacción entre los factores, a pesar de que el gráfico de interacción parece indicar que dicho efecto es irrelevante.
¿Qué efectos son relevantes para explicar el contenido promedio de grasa (en porcentaje) en la leche del ganado? Justifica tu respuesta.
¿Cuál es el modelo estimado?
Butterfat ~ Bread.
(Fc(4,95) = 49.8, p < .0001) ¿Consideras que se ha finalizado con el proceso de estimación del modelo? se rechaza Ho. Existe evidencia con un α = .05 que indica que las medias de Butterfat en Bread son diferentes. ¿Qué otra información deberíamos extraer?
# Ajuste del modelo
ajuste <- lm(Butterfat ~ Bread*Age,
data = ejer09)
# Selección del modelo
ajuste <- step(ajuste)## Start: AIC=-165.92
## Butterfat ~ Bread * Age
##
## Df Sum of Sq RSS AIC
## - Bread:Age 4 0.51387 16.094 -170.67
## <none> 15.580 -165.92
##
## Step: AIC=-170.67
## Butterfat ~ Bread + Age
##
## Df Sum of Sq RSS AIC
## - Age 1 0.274 16.368 -170.987
## <none> 16.094 -170.672
## - Bread 4 34.321 50.415 -64.487
##
## Step: AIC=-170.99
## Butterfat ~ Bread
##
## Df Sum of Sq RSS AIC
## <none> 16.368 -170.987
## - Bread 4 34.321 50.689 -65.946# Tabla ANOVA
anova(ajuste)# Análsis inferencial sobre los coeficientes del modelo
tidy(ajuste)Realizamos le test de comparaciones múltiples.
¿Qué conclusiones podemos extraer del análisis realizado? que las medias de tratamientos osn diferentes. ¿Cómo podrías representar la información obtenida para una mejor visualización de los resultados?
podemos calacular su homogeneidad de varianza , su modelo de ajuste de normalidad.
# Comparaciones múltiples
TukeyHSD(aov(Butterfat ~ Bread,
data = ejer09))## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Butterfat ~ Bread, data = ejer09)
##
## $Bread
## diff lwr upr p adj
## Canadian-Ayrshire 0.3785 0.01348598 0.74351402 0.0381804
## Guernsey-Ayrshire 0.8900 0.52498598 1.25501402 0.0000000
## Holstein-Fresian-Ayrshire -0.3905 -0.75551402 -0.02548598 0.0297910
## Jersey-Ayrshire 1.2325 0.86748598 1.59751402 0.0000000
## Guernsey-Canadian 0.5115 0.14648598 0.87651402 0.0016669
## Holstein-Fresian-Canadian -0.7690 -1.13401402 -0.40398598 0.0000007
## Jersey-Canadian 0.8540 0.48898598 1.21901402 0.0000000
## Holstein-Fresian-Guernsey -1.2805 -1.64551402 -0.91548598 0.0000000
## Jersey-Guernsey 0.3425 -0.02251402 0.70751402 0.0767002
## Jersey-Holstein-Fresian 1.6230 1.25798598 1.98801402 0.0000000-Diagnóstico del modelo Comenzamos con el diagnóstico del modelo ajustado.
¿Qué conclusiones extraemos del diagnóstico realizado? no hay homogeneidad de varianza Bartlett test of homogeneity of variances=0.0004, hay normalidad., puesto que no se rechaza ho.No existe evidencia con un α=.05 que indique que las varianzas de Burtterfat en Bread no sean diferentes. ¿Cómo propones que deberíamos modificar nuestro modelo? podemos hacer la transformacion de Boxcox o tranformacion Logaritmo
# Diagnóstico
fitinf <- fortify(ajuste)
# Varianza constante
bartlett.test(Butterfat ~ Bread, data = ejer09)##
## Bartlett test of homogeneity of variances
##
## data: Butterfat by Bread
## Bartlett's K-squared = 20.251, df = 4, p-value = 0.0004455# Normalidad (Obtenemos los pvalores del test de normalidad)
shapiro.test(fitinf$.stdresid[fitinf$Bread == "Ayrshire"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Bread == "Ayrshire"]
## W = 0.95688, p-value = 0.4836shapiro.test(fitinf$.stdresid[fitinf$Bread == "Canadian"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Bread == "Canadian"]
## W = 0.94549, p-value = 0.3037shapiro.test(fitinf$.stdresid[fitinf$Bread == "Guernsey"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Bread == "Guernsey"]
## W = 0.97024, p-value = 0.76shapiro.test(fitinf$.stdresid[fitinf$Bread == "Holstein-Fresian"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Bread == "Holstein-Fresian"]
## W = 0.90241, p-value = 0.04577shapiro.test(fitinf$.stdresid[fitinf$Bread == "Jersey"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Bread == "Jersey"]
## W = 0.95575, p-value = 0.4627¿Qué transformación debemos realizar sobre la variable respuesta? Transformamos la variable y ajustamos de nuevo el modelo completo. transformacion Logaritma en variable respuesta y variable independiente
# Transformación de boxcox
boxcox(ajuste) 
¿Qué modelo resulta tras la transformación de la respuesta? Butterfatinv ~ Bread + Age, donde Butterfatinv = 1/Butterfat (Fc(4,94) = 62.8, p < .0001) se rechaza Ho. Existe evidencia con un α = .05 que indica que las medias de Butterfat en Bread son diferentes.
¿Cuáles son las ecuaciones del modelo ajustado? Butterfatinv ~ 0.24-0.02BreadCanadian-0.04BreadGuernsey+0.02BreadHolstein-Fresian-0.05BreadJersey-0.00AgeMature
¿Cómo interpretamos los coeficientes de este modelo? Existe evidencia con un α = .05 que indica que las medias de Butterfat en Bread son diferentes A la vista de la tabla ANOVA ¿cuál debería ser el modelo a utilizar? Butterfatinv ~ Bread + Age
# Creamos la variable transformada
ejer09 <- ejer09 %>% mutate(Butterfatinv = 1/Butterfat)
#####################################
# Realizamos el ajuste de nuevo
ajuste <- lm(Butterfatinv ~ Bread*Age,
data = ejer09)
# Selección del modelo
ajuste <- step(ajuste)## Start: AIC=-784.54
## Butterfatinv ~ Bread * Age
##
## Df Sum of Sq RSS AIC
## - Bread:Age 4 0.0010253 0.033084 -789.39
## <none> 0.032059 -784.54
##
## Step: AIC=-789.39
## Butterfatinv ~ Bread + Age
##
## Df Sum of Sq RSS AIC
## <none> 0.033084 -789.39
## - Age 1 0.000698 0.033782 -789.30
## - Bread 4 0.087973 0.121057 -667.67# Tabla ANOVA
anova(ajuste)# Análsis inferencial sobre los coeficientes del modelo
tidy(ajuste)Realizamos el diagnóstico para este nuevo modelo.
¿Es adecuado el modelo obtenido? Justifica tu respuesta.si
se rechaza Ho. Existe evidencia con un α = .05 que indica que las medias de Butterfat en Bread son diferentes. El modelo ajustado para butterfatunc se cumple los supuestos de normalidad hapiro-Wilk normality test=0.628,homogeneidad de varianza Bartlett=0.80.
# Realizamos el ajuste de nuevo
ajuste <- lm(Butterfatinv ~ Bread, data = ejer09)
# Diagnóstico
fitinf <- fortify(ajuste)
# Varianza constante
bartlett.test(Butterfatinv ~ Bread, data = ejer09)##
## Bartlett test of homogeneity of variances
##
## data: Butterfatinv by Bread
## Bartlett's K-squared = 1.6238, df = 4, p-value = 0.8045# Normalidad (Obtenemos los pvalores del test de normalidad)
shapiro.test(fitinf$.stdresid[fitinf$Bread == "Ayrshire"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Bread == "Ayrshire"]
## W = 0.94642, p-value = 0.316shapiro.test(fitinf$.stdresid[fitinf$Bread == "Canadian"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Bread == "Canadian"]
## W = 0.96411, p-value = 0.6288shapiro.test(fitinf$.stdresid[fitinf$Bread == "Guernsey"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Bread == "Guernsey"]
## W = 0.95323, p-value = 0.4187shapiro.test(fitinf$.stdresid[fitinf$Bread == "Holstein-Fresian"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Bread == "Holstein-Fresian"]
## W = 0.94554, p-value = 0.3044shapiro.test(fitinf$.stdresid[fitinf$Bread == "Jersey"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Bread == "Jersey"]
## W = 0.97145, p-value = 0.7851-Predicción Construimos los intervalos de confianza para la predicción correspondiente al modelo ajustado.
¿Qué podemos concluir a la vista del gráfico de predicción? que las medias son diferentes.
¿Sería posible clasificar una oveja sin mirarla conociendo el contenido promedio de grasa de su leche? no. Justifica como procederías para contestar a esta cuestión.
# Predicción
# Secuencia de valores de predicción
newdata <- data.frame(Bread = c("Ayrshire","Canadian","Guernsey"
,"Holstein-Fresian","Jersey"))
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata,
interval="confidence"))
# Calculamos valore de predicción en la escala original
newdata$fit <- 1/newdata$fit
newdata$lwr <- 1/newdata$lwr
newdata$upr <- 1/newdata$upr
#Gráfico de predicción
ggplot(newdata, aes(x = fct_reorder(Bread,fit), y= fit)) +
geom_errorbar(aes(ymin=lwr, ymax=upr), width=.1) +
geom_line() +
geom_point() +
coord_flip() +
labs(y = "Contenido promedio de grasa", x = "Raza") +
theme_bw()## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
-\(Caso 9 (ANCOVA)\)
Los datos que se presentan son los pesos al nacer (en gramos) y las edades gestacionales estimadas (en semanas) de 12 bebés hombres y mujeres nacidos en un determinado hospital. Las edades medias son casi las mismas para ambos sexos. Se está interesado en estimar el peso de los bebes a partir de su sexo y su edad gestacional. Tenemos una variable predictora de tipo numérico (edad gestacional) y otra de tipo categórico (sexo). Consideramos un modelo ANCOVA con un regresor y un factor. Como el sexo viene codificado como 1 (Boy) y 2 (Girl), construimos el factor asociado antes de la representación gráfica de los datos.
# Lectura de datos
previo <- read_csv("https://goo.gl/B3yoLJ", col_types = "ddc")
# Recodificación del factor
ejer10 <- previo %>% mutate(sex=fct_recode(sex,"Boy"="1","Girl"="2"))Para la representación gráfica de los datos consideramos las tres posibilidades de modelización para estos datos:
Una única recta de regresión entre peso y edad. Se considera irrelevante la asociación entre el sexo y edad, y el posible efecto del sexo del bebé en su peso final. Dos rectas de regresión paralelas (una por cada sexo). Se considera irrelevante la asociación entre el sexo y edad, pero se considera un posible efecto del sexo del bebé en su peso final. Dos rectas distintas (una por cada sexo). Se considera relevante la asociación entre el sexo y edad para el peso final del bebé.
¿Hay evidencias de cuál puede ser el mejor modelo de los tres planteados? los tres graficos casi tiene comportamiento similares,entonces se puede hacer el modelo ,para sabeer si la variable es explicada a la variable respuesta. Justifica tu respuesta.
# Representación gráfica de los modelos de regresión
# Modelo con una única recta
M0 <- lm(weight ~ age,data = ejer10)
# M1: modelo con rectas paralelas
M1 <- lm(weight ~ sex + age,data = ejer10)
# M2: modelo con rectas no paralelas
M2 <- lm(weight ~ sex * age,data = ejer10)
# grid de valores para construir los modelos
grid <- ejer10 %>% data_grid(sex,age) %>%
gather_predictions(M0,M1,M2)
# Gráfico
ggplot(ejer10,aes(age,weight,colour = sex)) +
geom_point() +
geom_line(data = grid, aes(y = pred)) +
facet_wrap(~ model) +
labs(x = "Edad Gestacional", y = "Peso") +
theme_bw()
-Estimación del modelo Comparamos los modelos construidos en el análisis preliminar para decidirnos entre uno de ellos. Para ello partimos del modelo más complejo y realizamos el proceso de selección de efectos.
¿Qué efectos están presentes en el modelo según el proceso de selección automático? Si comparamos este resultado con la tabla ANOVA resultante Existe evidencia con un α = .05 que indica que las medias weight en aged son diferentes. Existe evidencia con un α = .05 que indica que las medias weight en no son diferentes.
¿el modelo seleccionado sería el mismo? ¿Qué opción te parece más adecuada? ¿qué implicaciones experimentales tienen ambos modelos? Si consideramos como válido el proporcionado por el procedimiento automático de selección ¿Cómo interpretamos los coeficientes de dicho modelo? ¿Podríamos eliminar la interceptación de este modelo?
# Selección del modelo
ajuste <- M2
ajuste <- step(ajuste)## Start: AIC=253.05
## weight ~ sex * age
##
## Df Sum of Sq RSS AIC
## - sex:age 1 6346.2 658771 251.28
## <none> 652425 253.05
##
## Step: AIC=251.28
## weight ~ sex + age
##
## Df Sum of Sq RSS AIC
## <none> 658771 251.28
## - sex 1 157304 816074 254.42
## - age 1 1094940 1753711 272.78# Tabla ANOVA
anova(ajuste)# Coeficientes del modelo
tidy(ajuste)Diagnóstico del modelo Realizamos el diagnóstico del modelo que contiene los efectos principales de edad y sexo del bebé. Realizamos tanto los gráficos de residuos como los tests de diagnóstico.
¿Se detecta algún problema con las hipótesis del modelo? ¿Cómo afecta eso al modelo construido?
# Diagnóstico del modelo
fitinf <- fortify(ajuste)
# Residuos versus ajustados
ggplot(fitinf, aes(x = .fitted, y =.stdresid)) +
geom_point() +
facet_wrap(~ sex, scales="free_x") +
theme_bw()
# Residuos versus predictora
ggplot(fitinf, aes(x = age, y =.stdresid)) +
geom_point() +
facet_wrap(~ sex, scales="free_x") +
theme_bw()
# Tests para verificar las hipótesis
# Varianza constante
bartlett.test(.stdresid ~ sex, data = fitinf)##
## Bartlett test of homogeneity of variances
##
## data: .stdresid by sex
## Bartlett's K-squared = 0.64847, df = 1, p-value = 0.4207# Normalidad (Obtenemos los pvalores del test de normalidad)
shapiro.test(fitinf$.stdresid[fitinf$sex == "Boy"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$sex == "Boy"]
## W = 0.93081, p-value = 0.3888shapiro.test(fitinf$.stdresid[fitinf$sex == "Girl"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$sex == "Girl"]
## W = 0.87317, p-value = 0.0717-Predicción Creamos un grid de valores para la edad gestacional y construimos los modelos de predicción asociados con cada nivel del factor.
¿Qué conclusiones podemos extraer de los modelos de predicción construidos? ¿Qué diferencia de peso estimamos para una edad gestacional de 38 semanas entre ambos sexos? ¿Podemos tomar este modelo de predicción como válido para futuros embarazos? ¿Cómo procederíamos y que dificultades nos podríamos encontrar?
newdata <- data.frame(age = rep(seq(min(ejer10$age),
max(ejer10$age), .01),
each=2),
sex = factor(c("Boy", "Girl")))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(ajuste, newdata,
interval="confidence"))
# Gráfico de la predicción
ggplot(newdata, aes(x = age, y = fit, color = sex)) +
geom_line() +
geom_ribbon(aes(ymax = upr, ymin = lwr, fill = sex), alpha = 1/5) +
labs(x = "Edad gestacional", y = "Peso al nacer",
title = "Bandas de predicción") +
theme_bw()
-\(Caso 10 (ANCOVA)\) Se lleva a cabo una investigación sobre diversas malformaciones del sistema nervioso central registradas en nacidos vivos en Gales del Sur, Reino Unido. El estudio fue diseñado para determinar el efecto de la dureza del agua sobre la incidencia de tales malformaciones. La información registrada son: NoCNS = recuento de nacimientos sin problema CNS; An = conteo de nacimientos de Anencephalus; Sp = conteo de nacimientos de espina bífida; Otro = recuento de otros nacimientos del SNC; Agua = endurecimiento del agua; Trabajo = un factor con niveles Manual o No manual en función del tipo de trabajo realizado por los padres.
# Lectura de datos
previo <- read_csv("https://goo.gl/bNOSxt", col_types = "cdddddc")
# Calculamos el número total de malformaciones para utilizarla como variable
ejer11 <- previo %>% mutate(CNS=An+Sp+Other)Dado que estamos interesados en establecer un modelo para el número total de malformaciones, deberíamos construir en primer lugar una variable que contenga dichos valores como la suma de todos los tipos de malformaciones registradas. Esta nueva variable se utilizará como respuesta en el modelo propuesto. Como posibles variables predictoras tenemos el tipo de trabajo de los padres y la dureza del agua. Planteamos un modelo ANCOVA con un regresor y un factor. No consideramos el distrito origen de los sujetos para este modelo. Procedemos con el caso anterior para representar las diferentes opciones de modelización para estos datos. ¿Qué información proporcionan los gráficos obtenidos de cara a la elección del posible modelo final sobre estos datos?
# Representación gráfica de los modelos de regresión
# Modelo con una única recta
M0 <- lm(CNS ~ Water,data = ejer11)
# M1: modelo con rectas paralelas
M1 <- lm(CNS ~ Work + Water,data = ejer11)
# M2: modelo con rectas no paralelas
M2 <- lm(CNS ~ Work * Water,data = ejer11)
# grid de valores para construir los modelos
grid <- ejer11 %>% data_grid(Work,Water) %>%
gather_predictions(M0,M1,M2)
# Gráfico
ggplot(ejer11,aes(Water,CNS,colour = Work)) +
geom_point() +
geom_line(data = grid, aes(y = pred)) +
facet_wrap(~ model) +
labs(x = "Dureza del agua", y = "Número malformaciones") +
theme_bw()
-Estimación del modelo Seleccionamos el mejor modelo para estos datos partiendo del modelo más complejo (interacción entre dureza del agua y tipo de trabajo). ¿Qué tipo de modelo se ha estimado finalmente? ¿Cuales son las rectas de regresión estimadas? ¿Cómo interpretamos los coeficientes de dichos modelos? ¿Qué sentido tiene la interceptación del modelo?
# Selección del modelo
ajuste <- M2
ajuste <- step(ajuste)## Start: AIC=105.6
## CNS ~ Work * Water
##
## Df Sum of Sq RSS AIC
## <none> 7135.0 105.60
## - Work:Water 1 1848.8 8983.8 107.29# Tabla ANOVA
anova(ajuste)# Coeficientes del modelo
tidy(ajuste)Diagnóstico del modelo Realizamos el diagnóstico gráfico y basado en los tests de diagnóstico para este tipo de modelos.
¿Qué indica el gráfico de residuos versus ajustados para cada nivel del factor? ¿Cómo interpretamos los resultados de los tests de diagnóstico? ¿Cómo procedemos a partir de ahora?
# Diagnóstico del modelo
fitinf <- fortify(ajuste)
# Residuos versus ajustados
ggplot(fitinf, aes(x = .fitted, y =.stdresid)) +
geom_point() +
facet_wrap(~ Work, scales="free_x") +
theme_bw()
# Residuos versus predictora
ggplot(fitinf, aes(x = Water, y =.stdresid)) +
geom_point() +
facet_wrap(~ Work, scales="free_x") +
theme_bw()
# Tests para verificar las hipótesis
# Varianza constante
bartlett.test(.stdresid ~ Work, data = fitinf)##
## Bartlett test of homogeneity of variances
##
## data: .stdresid by Work
## Bartlett's K-squared = 11.677, df = 1, p-value = 0.0006327# Normalidad (Obtenemos los pvalores del test de normalidad)
shapiro.test(fitinf$.stdresid[fitinf$Work == "Manual"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Work == "Manual"]
## W = 0.92992, p-value = 0.5154shapiro.test(fitinf$.stdresid[fitinf$Work == "NonManual"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Work == "NonManual"]
## W = 0.93462, p-value = 0.559Probamos con las transformaciones de Box-Cox.
¿Existe una transformación adecuada para este modelo? Transformamos la respuesta y reajustamos el modelo.
boxcox(ajuste)
¿Qué modelo resulta del nuevo proceso de estimación? ¿Cómo interpretamos los coeficientes de dicho modelo?
# Creamos la variable transformada
ejer11 <- ejer11 %>% mutate(CNSlog = log(CNS))
#################################################
# Nuevo ajuste
ajuste <- lm(CNSlog ~ Work * Water,data = ejer11)
ajuste <- step(ajuste)## Start: AIC=-16.94
## CNSlog ~ Work * Water
##
## Df Sum of Sq RSS AIC
## - Work:Water 1 0.25147 3.6176 -17.788
## <none> 3.3661 -16.941
##
## Step: AIC=-17.79
## CNSlog ~ Work + Water
##
## Df Sum of Sq RSS AIC
## - Water 1 0.4770 4.0946 -17.8066
## <none> 3.6176 -17.7884
## - Work 1 7.1815 10.7991 -2.2901
##
## Step: AIC=-17.81
## CNSlog ~ Work
##
## Df Sum of Sq RSS AIC
## <none> 4.0946 -17.8066
## - Work 1 7.1815 11.2761 -3.5985# Tabla ANOVA
anova(ajuste)# Coeficientes del modelo
tidy(ajuste)¿Qué ocurre con las hipótesis de este modelo?
# Nuevo diagnóstico
fitinf <- fortify(ajuste)
# Varianza constante
bartlett.test(.stdresid ~ Work, data = fitinf)##
## Bartlett test of homogeneity of variances
##
## data: .stdresid by Work
## Bartlett's K-squared = 0.92472, df = 1, p-value = 0.3362# Normalidad (Obtenemos los pvalores del test de normalidad)
shapiro.test(fitinf$.stdresid[fitinf$Work == "Manual"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Work == "Manual"]
## W = 0.89776, p-value = 0.2758shapiro.test(fitinf$.stdresid[fitinf$Work == "NonManual"])##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid[fitinf$Work == "NonManual"]
## W = 0.92223, p-value = 0.4483Predicción Construimos los intervalos de confianza para la predicción correspondiente a cada uno de los niveles del factor. Por el momento seguimos trabajando con la variable transformada.
¿Cómo obtendríamos los valores de predicción en la escala original del número de malformaciones?
# Secuencia de valores de predicción
newdata <- data.frame(Work = unique(ejer11$Work))
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata,
interval="confidence"))
# Gráfico resultante (predicción y bandas de predicción).
# Ordenamos las clases de acuerdo al valor de fit
ggplot(newdata, aes(x = fct_reorder(Work,fit), y= fit)) +
geom_errorbar(aes(ymin=lwr, ymax=upr), width=.1) +
geom_line() +
geom_point() +
labs(x = "Tipo de trabajo", y = "Logaritmo número de malformaciones") +
coord_flip() +
theme_bw()## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
-\(Caso 11 (ANCOVA)\) Se ha realizado un estudio para establecer la calidad de los vinos de la variedad Pino Noir en función de un conjunto de características analizadas. Las características analizadas son claridad, aroma, cuerpo, olor y matiz. Para medir la calidad se organiza una cata ciega a un conjunto de expertos y se calcula la puntuación final de cada vino a partir de la información de todos ellos. Además se registra la región de procedencia del vino por si puede influir en la calidad del vino.
# Lectura de datos
ejer12 <- read_csv("https://goo.gl/OX9wgM", col_types = "ddddddc")En este caso tenemos un modelo ANCOVA con un factor (región de procedencia) y cinco predictoras numéricas, para tratar de explicar la calidad del vino. En este caso no realizamos el análisis gráfico y pasamos directamente a la estimación del modelo.
¿Qué tipo de gráfico podríamos construir para representar la información de este modelo? ¿Nos aportaría información a la hora de plantear el posible modelo estadístico?
-Estimación del modelo Planteamos el modelo completo y realizamos el correspondiente proceso de selección automática de efectos.
¿Qué modelo resulta del proceso de estimación? ¿Cómo podemos interpretar los diferentes efecto de interacción que aparecen en él?
# Gráfico del modelo díficil de realizar
# Pasamos al ajuste directamente
ajuste <- lm(calidad ~ region * (claridad + aroma + cuerpo + olor + matiz),
data = ejer12)
ajuste <- step(ajuste)## Start: AIC=-0.12
## calidad ~ region * (claridad + aroma + cuerpo + olor + matiz)
##
## Df Sum of Sq RSS AIC
## - region:aroma 2 0.0229 14.712 -4.0591
## - region:olor 2 1.0100 15.699 -1.5913
## - region:matiz 2 1.5643 16.253 -0.2728
## <none> 14.689 -0.1184
## - region:cuerpo 2 4.1828 18.872 5.4032
## - region:claridad 2 4.7222 19.411 6.4741
##
## Step: AIC=-4.06
## calidad ~ region + claridad + aroma + cuerpo + olor + matiz +
## region:claridad + region:cuerpo + region:olor + region:matiz
##
## Df Sum of Sq RSS AIC
## - aroma 1 0.1698 14.882 -5.6230
## - region:olor 2 1.0633 15.775 -5.4073
## <none> 14.712 -4.0591
## - region:matiz 2 1.9162 16.628 -3.4065
## - region:cuerpo 2 4.4792 19.191 2.0408
## - region:claridad 2 4.7625 19.474 2.5978
##
## Step: AIC=-5.62
## calidad ~ region + claridad + cuerpo + olor + matiz + region:claridad +
## region:cuerpo + region:olor + region:matiz
##
## Df Sum of Sq RSS AIC
## - region:olor 2 1.0526 15.934 -7.0261
## <none> 14.882 -5.6230
## - region:matiz 2 2.4006 17.282 -3.9402
## - region:cuerpo 2 4.5777 19.459 0.5684
## - region:claridad 2 5.2764 20.158 1.9088
##
## Step: AIC=-7.03
## calidad ~ region + claridad + cuerpo + olor + matiz + region:claridad +
## region:cuerpo + region:matiz
##
## Df Sum of Sq RSS AIC
## <none> 15.934 -7.0261
## - region:matiz 2 3.0051 18.939 -4.4608
## - region:claridad 2 4.2682 20.203 -2.0074
## - region:cuerpo 2 4.3696 20.304 -1.8172
## - olor 1 9.7000 25.634 9.0412tidy(ajuste)-Diagmóstico del modelo ¿Cómo deberíamos proceder para realizar el diagnóstico de este modelo?
-Predicción En este caso el gráfico de predicción resulta muy complejo debido a la gran cantidad de efecto presentes en el modelo. Nos podemos limitar a obtener una tabla de predicciones para valores específicos de las predictoras.
¿Qué valores de las predictoras propones para este proceso? ¿Cómo realizaríamos este proceso? La opción más adecuada en esta situación es el establecimiento de un análisis de escenarios para resolver este problema. En un análisis de escenarios propone tres valores de la variable numérica (uno bajo, otro intermedio y otro alto) y predice los valores correspondientes de la respuesta. ¿Como actuaríamos en este problema?
-\(Caso 12 (ANCOVA)\) El banco de datos de Puromycin contiene 23 mediciones sobre la velocidad de reacción enzimática frente a la concentración de sustrato para células tratadas o no tratadas con Puromicina. Las variables registradas son:
conc: Concentración de sustrato en partes por millón (ppm). rate: Velocidad instántanea de reacción (recuentos/min/min). state: Estado (Tratatado o no tratado con Puromicina)
# Carga de datos de librería de R
data("Puromycin")Pretendemos construir un modelo que nos permita conocer la velocidad de reacción en función de la concentración del sustrato y de si ha sido tratado o no con Puromicina. Se trata de un modelo ANCOVA con una factor y una variable numérica. Comenzamos con la representación gráfica del modelo de interacción. ¿Consideras adecuado un modelo lineal para este diseño experimnetal? ¿Qué opciones de modelización podríamos plantear?
# Gráfico
ggplot(Puromycin) +
geom_point(aes(conc,rate, colour = state)) +
geom_smooth(aes(conc,rate, colour = state),
method= mgcv::gam, se=FALSE) +
labs(x = "Concengración", y = "Velocidad de reacción") +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
Realizamos un nuevo gráfico con las ecuaciones de suavizado
# Gráfico
ggplot(Puromycin) +
geom_point(aes(conc,rate, colour = state)) +
geom_smooth(aes(conc,rate, colour = state),
method= loess, se=FALSE) +
labs(x = "Concentración", y = "Velocidad de reacción") +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
Estimación del modelo Construimos el modelo de suavizado para las predictoras consideradas.
¿Cuáles son las principales consecuencias del modelo ajustado? ¿Como interpretamos el coeficiente asociado con si el sujeto ha sido o no tratado? ¿Consideras necesaria la eliminación de algún efecto en el modelo porpuesto? ¿Podemos justificar asumir este modelo frente a otro donde no consideramos la interacción entre el estado y la concentración? Justifica tu respuesta.
ajuste <- gam(rate ~ state + s(conc, k = 10, m = 2, bs = "ps", by = state),
data = Puromycin)## Warning in smooth.construct.ps.smooth.spec(object, dk$data, dk$knots): basis
## dimension is larger than number of unique covariatessummary(ajuste)##
## Family: gaussian
## Link function: identity
##
## Formula:
## rate ~ state + s(conc, k = 10, m = 2, bs = "ps", by = state)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 138.891 2.711 51.227 < 2e-16 ***
## stateuntreated -26.020 3.930 -6.621 1.37e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(conc):statetreated 3.913 4.151 82.75 < 2e-16 ***
## s(conc):stateuntreated 3.563 3.858 37.45 1.08e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.961 Deviance explained = 97.6%
## GCV = 148.62 Scale est. = 87.386 n = 23Ajustamos ambos modelos y los comparamos. ¿Qué modelo debemos asumir como definitivo?
ajuste <- gam(rate ~ state + s(conc, k = 10, m = 2, bs = "ps", by = state),
data = Puromycin)## Warning in smooth.construct.ps.smooth.spec(object, dk$data, dk$knots): basis
## dimension is larger than number of unique covariatesajuste2 <- gam(rate ~ state + s(conc, k = 10, m = 2, bs = "ps"),
data = Puromycin)## Warning in smooth.construct.ps.smooth.spec(object, dk$data, dk$knots): basis
## dimension is larger than number of unique covariatesdata.frame(ajuste$gcv.ubre,ajuste2$gcv.ubre)-Diagnóstico del modelo Procedemos con el diagnóstico del modelo. Comenzamos con el análisis gráfico y de influencia.
¿Qué podemos concluir del análisis realizado?
gam.check(ajuste)
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 8 iterations.
## The RMS GCV score gradient at convergence was 2.925858e-06 .
## The Hessian was positive definite.
## Model rank = 20 / 20
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(conc):statetreated 9.00 3.91 1.52 0.99
## s(conc):stateuntreated 9.00 3.56 1.52 0.99-Predicción Obtenemos el gráfico de predicción asociado con el modelo ajustado.
¿Cómo valoras la predicción obtneida? Si el modelo cumple con las hipótesis ¿por qué la predicción parece funcionar tan mal? ¿qué podemos hacer ahora?
# Creamos grid de predicción
# Utlizamos la opción each para indicar el número de repeticiones (tantas como niveles del factor)
newdata <- data.frame(conc = rep(seq(min(Puromycin$conc), max(Puromycin$conc), .01),each=2), state = factor(c("treated", "untreated")))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(ajuste, newdata, type = "response", se.fit=TRUE))
# Calculamos intervalos de confianza
newdata <- newdata %>% mutate(
upr = fit + (1.96 * se.fit),
lwr = fit - (1.96 * se.fit)
)
# Gráfico de la predicción
ggplot(newdata, aes(x = conc, y = fit, color = state)) +
geom_line() +
geom_ribbon(aes(ymax = upr, ymin = lwr, fill = state), alpha = 1/5) +
labs(x = "Concentración", y = "Velocidad de reacción", title = "Bandas de predicción") +
theme_bw()
¿Qué tipo de modelo hemos ajusatdo? ¿Cómo valoras la capacidad explicativa de este modelo? Justifica tu respuesta
ajuste <- lm(rate ~ state * (conc + I(conc^2)),
data = Puromycin)
ajuste <- step(ajuste)## Start: AIC=132.91
## rate ~ state * (conc + I(conc^2))
##
## Df Sum of Sq RSS AIC
## - state:I(conc^2) 1 278.82 4691.9 132.32
## <none> 4413.1 132.91
## - state:conc 1 581.45 4994.6 133.75
##
## Step: AIC=132.32
## rate ~ state + conc + I(conc^2) + state:conc
##
## Df Sum of Sq RSS AIC
## <none> 4691.9 132.32
## - state:conc 1 882.8 5574.7 134.28
## - I(conc^2) 1 9210.1 13902.0 155.30tidy(ajuste)glance(ajuste)¿Cómo interpretamos los resultados de los tests de diagnóstico? ¿Cómo procedemos a partir de ahora? ¿Dispones de modelos adecuados para este tipo de datos?
# Varianza constante
bptest(ajuste)##
## studentized Breusch-Pagan test
##
## data: ajuste
## BP = 3.3308, df = 4, p-value = 0.5041# Normalidad
shapiro.test(fitinf$.stdresid)##
## Shapiro-Wilk normality test
##
## data: fitinf$.stdresid
## W = 0.93588, p-value = 0.3016# Independencia
dwtest(ajuste)##
## Durbin-Watson test
##
## data: ajuste
## DW = 1.2334, p-value = 0.001971
## alternative hypothesis: true autocorrelation is greater than 0¿Qué ocurre si transformamos tanto la concentración como velocidad mediante el logaritmo y tratamos de ajustar un modelo a esos datos?
-\(Caso 13 (ANCOVA)\)
El banco de datos presenta la información referida al nacimiento y mortalidad infantil de 800 niños nacidos en el estado de Carolina del Norte. Las variables consideradas en el estudio son:
plural: Número de hijos nacidos del embarazo. sex: Sexo del bebe. mage: Edad de la madre. weeks: Semanas completas de gestación. marital: Estado matrimonial (“married”=1; “not married”=2). racemom: Raza de la madre (“other non white”=0,“White”=1,“Black”=2,“America indian”=3,“Chinese”=4,“Hawaiian”=5,“Filipino”=6,“Other asian”=7). hispmom: Madre de origen hispánico (“Cuban”=C,“Mexican”=M,“Non-Hispanic”=N,“Other”=O,“Puerto Rican”=P,“Central/South american”=S,“Not classificable”=U). gained: Peso ganado durante el embarazo (en libras). smoke: Madre fumadora (“Yes”=1,“No”=0). drink: Madre bebedora (“Yes”=1,“No”=0). tounces: Peso del bebe (en onzas). tgrams: Peso del bebe (en gramos). low: Bebe de poco peso (“Yes”=1,“No”=0). premie: Bebe prematuro (“Yes”=1,“No”=0). Obtén un modelo que trate de explicar el peso del bebé al nacer (en gramos) en función del sexo del bebe, de la edad de la madre, de las semansa completas de gestación, de si la madre es bebedora o fumadora, y de si el bebé ha sido prematuro.
A continuación tienes el código para la lectura de datos y la codificación de los factores.
ncbirth <- read_csv("https://goo.gl/mB9Jcn", col_types = "dcddcccdccddcc")
str(ncbirth)## spc_tbl_ [800 × 14] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ plural : num [1:800] 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : chr [1:800] "1" "2" "1" "1" ...
## $ mage : num [1:800] 32 32 27 27 25 28 25 15 37 21 ...
## $ weeks : num [1:800] 40 37 39 39 39 43 39 42 41 39 ...
## $ marital: chr [1:800] "1" "1" "1" "1" ...
## $ racemom: chr [1:800] "1" "1" "1" "1" ...
## $ hispmom: chr [1:800] "N" "N" "N" "N" ...
## $ gained : num [1:800] 38 34 12 15 32 32 75 25 31 28 ...
## $ smoke : chr [1:800] "0" "0" "0" "0" ...
## $ drink : chr [1:800] "0" "0" "0" "0" ...
## $ tounces: num [1:800] 111 116 138 136 121 117 143 113 139 120 ...
## $ tgrams : num [1:800] 3147 3289 3912 3856 3430 ...
## $ low : chr [1:800] "0" "0" "0" "0" ...
## $ premie : chr [1:800] "0" "0" "0" "0" ...
## - attr(*, "spec")=
## .. cols(
## .. plural = col_double(),
## .. sex = col_character(),
## .. mage = col_double(),
## .. weeks = col_double(),
## .. marital = col_character(),
## .. racemom = col_character(),
## .. hispmom = col_character(),
## .. gained = col_double(),
## .. smoke = col_character(),
## .. drink = col_character(),
## .. tounces = col_double(),
## .. tgrams = col_double(),
## .. low = col_character(),
## .. premie = col_character()
## .. )
## - attr(*, "problems")=<externalptr># Recodificación del factor
ncbirth <- ncbirth %>%
mutate(sex=fct_recode(sex,"male"="1","female"="2"),
marital=fct_recode(marital,"married"="1","not married"="2"),
racemom=fct_recode(racemom,"other non white"="0",
"White"="1",
"Black"="2",
"America indian"="3",
"Chinese"="4",
"Hawaiian"="5",
"Filipino"="6",
"Other asian"="7",
"Other"="8"),
hispmom=fct_recode(hispmom,"Cuban"="C",
"Mexican"="M",
"Non-Hispanic"="N",
"Other"="O",
"Puerto Rican"="P",
"Central/South american"="S",
"U"="Not classificable"),
smoke=fct_recode(smoke,"Yes"="1","No"="0"),
drink=fct_recode(drink,"Yes"="1","No"="0"),
low=fct_recode(low,"Yes"="1","No"="0"),
premie=fct_recode(premie,"Yes"="1","No"="0"))## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `racemom = fct_recode(...)`.
## Caused by warning:
## ! Unknown levels in `f`: 0, 5, 6
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 1 remaining warning.-\(Caso 14 (ANOVA)\) Un grupo de varones adultos con edades comprendidas entre 30 y 65 años participaron en un estudio para investigar la relación entre el consumo de carne y el colesterol. Los sujetos fueron organizados en tres grupos de acuerdo a tres dietas diferentes con una duración de 20 semanas. Las variables consideradas son:
SUBJ: Sujeto. Dieta: Tipo de dieta (“BEEF” = carne de vaca unicamente, “PORK” = carne de cerdo unicamente, “CHFISH” = carne de pollo y pescado unicamente). chol: Nivel de colesterol. Obtén un modelo que trate de explicar el nivel de colesterol en función de la dieta seguida.
A continuación tienes el código para la lectura de datos.
serumcho <- read_csv("https://goo.gl/ghxka2", col_types = "iddd")
str(serumcho)## spc_tbl_ [347 × 4] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ SUBJ : int [1:347] 1 2 3 4 5 6 7 8 9 10 ...
## $ BEEF : num [1:347] 241 218 261 190 238 256 248 224 225 238 ...
## $ PORK : num [1:347] 245 197 199 162 191 182 160 180 208 227 ...
## $ CHFISH: num [1:347] 249 222 221 215 207 193 205 227 203 180 ...
## - attr(*, "spec")=
## .. cols(
## .. SUBJ = col_integer(),
## .. BEEF = col_double(),
## .. PORK = col_double(),
## .. CHFISH = col_double()
## .. )
## - attr(*, "problems")=<externalptr># Reorganizamos los datos
caso14 <- gather(serumcho,`BEEF`,`PORK`,`CHFISH`, key = "dieta", value = colesterol)