Taller IPERT - R Avanzado
D. S. Fernandez del Viso
1/31/2020
Taller R Avanzado: Regresión Múltiple, ANOVA, GLM
Lista de paquetes a instalar
- readxl
- car
- EnvStats
- tidyverse
- coefplot
- leaps
- rmarkdown
- ggplot2
- MASS
- multcomp
- effects
- gplots
- agricolae
- plyr
- caret
Lista de archivos a usar
Los siguientes archivos deben estar disponible en su directorio de trabajo de R:
Objetivos Generales
Conocer y aplicar los métodos disponibles en R para obtener los parámetros de modelos de regresión múltiple, evaluar su idoneidad para describir los datos, e interpretar su significado.
Aplicar funciones de R para diseños de investigación apropiados para el análisis de varianza.
Conocer y aplicar funciones de R para extender los modelos lineales a situaciones en las que la variable dependiente no tiene distribución normal (GLM).
DIA 1
Regresión Lineal Múltiple
Cuando tenemos más de una variable predictora (“independiente”), la regresión lineal simple viene a ser una regresión múltiple.
\[Y_j=\alpha+\beta_1X_{1j}+\beta_2X_{2j}+...+\beta_nX_{nj}\]
Ejemplo 1 y Datos
Usaremos como ejemplo los datos bmi en el archivo mod_empiricos.xlsx. Estos son datos de individuos adultos entre 21 y 79 años, con las siguientes variables: BMI, índice de masa corporal (\(kg/m^{2}\)); Age, edad (años); Cholesterol, niveles de colesterol en sangre, (\(mg/dL\)); Glucose, niveles de glucosa en la sangre, (\(mg/dL\)).
Para leer los datos emplearemos el siguiente código (obtenido del menú Import Dataset/From Excel):
library(readxl)
reg.multiple <- read_excel("Data/mod_empiricos.xlsx", sheet = "bmi")
# para ver las primeras seis filas de datos
head(reg.multiple)## # A tibble: 6 × 4
## BMI Age Cholesterol Glucose
## <dbl> <dbl> <dbl> <dbl>
## 1 19.3 21 178 95
## 2 24.5 57 250 98
## 3 24.7 46 176 102
## 4 47.9 47 171 105
## 5 44.2 61 222 101
## 6 29.9 74 156 72Pruebas de supuestos para regresión paramétrica usando el método de mínimos cuadrados (OLS, ordinary least square)
- Normalidad: Gráficamente: Q-Q plot; Estadísticamente: prueba de Shapiro & Wilk
- Independencia: los valores de \(Y_j\) son independientes entre si, lo asumimos (no hay autocorrelación).
- Linealidad: relación lineal entre variable dependiente y cada una de las independientes; gráficas individuales. Puede arreglarse con transformación; prueba de Box Tidwell.
- Homocedasticidad: varianza de la variable dependiente (residuales) no varía con los valores de las variables independientes; gráfica de los residuales.
- Multicolinealidad: las variables independientes no deben estar correlacionadas entre si.
Normalidad
- Evaluación gráfica de normalidad de la variable dependiente.
library(EnvStats)
EnvStats::qqPlot(reg.multiple$BMI, add.line = TRUE, points.col = "blue", line.col = "red")
- Prueba estadística de normalidad Shapiro & Wilk (\(H_0:distribución\enspace normal\))
shapiro.test(reg.multiple$BMI)##
## Shapiro-Wilk normality test
##
## data: reg.multiple$BMI
## W = 0.96315, p-value = 0.07547Autocorrelación
La prueba de Durbin-Watson nos permite evaluar si ocurre autocorrelación entre los valores residuales de la variable dependiente. Por ejemplo, una variable que depende del tiempo, presenta autocorrelación. La \(H_0\) es que la autocorrelación en los residuales del modelo es 0.
- Prueba de Durbin-Watson para autocorrelación
library(car)
model <- lm(BMI ~ Age + Cholesterol + Glucose, data = reg.multiple)
dbt <- durbinWatsonTest(model, simulate = TRUE)
dbt## lag Autocorrelation D-W Statistic p-value
## 1 0.06609505 1.823275 0.55
## Alternative hypothesis: rho != 0Linealidad
- Prueba de Box Tidwell para linealidad
library(car)
boxTidwell(BMI ~ Age + Cholesterol + Glucose, data = reg.multiple)## MLE of lambda Score Statistic (t) Pr(>|t|)
## Age -19.8117 -0.8984 0.3732
## Cholesterol 20.2797 0.1879 0.8517
## Glucose -2.3134 0.1057 0.9162
##
## iterations = 10
##
## Score test for null hypothesis that all lambdas = 1:
## F = 0.28618, df = 3 and 51, Pr(>F) = 0.8352Multicolinealidad
- Matriz de correlación y gráficas de puntos con líneas de regresión
library(car)
#matriz de correlación
cor(reg.multiple)## BMI Age Cholesterol Glucose
## BMI 1.0000000 0.1863333 0.28110580 0.22115998
## Age 0.1863333 1.0000000 0.16966704 0.26705577
## Cholesterol 0.2811058 0.1696670 1.00000000 0.08266412
## Glucose 0.2211600 0.2670558 0.08266412 1.00000000#gráfica de regresiones en parejas, con línea de regresión
scatterplotMatrix(reg.multiple, ~ BMI + Age + Cholesterol + Glucose,
smooth = list(lty = 2), id = TRUE,
regLine = list(lty = 1, col = "red"),
col = "blue")## Warning in applyDefaults(legend, defaults = list(coords = NULL, pt.cex = cex, :
## unnamed legend arguments, will be ignored
Modelos de regresión lineal múltiple
A continuación se calculan los parámetros de diversos modelos de regresión, y se incluye una prueba de homocedasticidad (homogeneidad de la varianza) para cada modelo.
Modelo completo (todas las variables)
modRM <- lm(BMI ~ Age + Cholesterol + Glucose,
data = reg.multiple)
summary(modRM)##
## Call:
## lm(formula = BMI ~ Age + Cholesterol + Glucose, data = reg.multiple)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6255 -4.5473 -0.8179 3.7439 18.8116
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.81510 5.07180 3.315 0.00164 **
## Age 0.04103 0.05631 0.729 0.46939
## Cholesterol 0.04819 0.02487 1.938 0.05791 .
## Glucose 0.01974 0.01493 1.322 0.19188
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.34 on 54 degrees of freedom
## Multiple R-squared: 0.127, Adjusted R-squared: 0.07855
## F-statistic: 2.62 on 3 and 54 DF, p-value: 0.06008aic <- AIC(modRM)
sprintf("AIC = %.2f", aic)## [1] "AIC = 384.68"# prueba de homocedasticidad (Non-constant Variance Score Test)
# Ho:la varianza es constante en el ámbito de la predicción de Y
ncvTest(modRM)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.02932531, Df = 1, p = 0.86403spreadLevelPlot(modRM)
##
## Suggested power transformation: 0.757141Modelo eliminando la variable Age, por ser la menos correlaciona con BMI, y tener la mayor correlación con Glucose
modRM <- lm(BMI ~ Cholesterol + Glucose,
data = reg.multiple)
summary(modRM)##
## Call:
## lm(formula = BMI ~ Cholesterol + Glucose, data = reg.multiple)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.125 -4.748 -1.134 3.321 18.832
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.95236 4.80502 3.736 0.000446 ***
## Cholesterol 0.05098 0.02447 2.083 0.041902 *
## Glucose 0.02254 0.01437 1.569 0.122449
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.313 on 55 degrees of freedom
## Multiple R-squared: 0.1185, Adjusted R-squared: 0.08641
## F-statistic: 3.696 on 2 and 55 DF, p-value: 0.0312aic <- AIC(modRM)
sprintf("AIC = %.2f", aic)## [1] "AIC = 383.25"# prueba de homocedasticidad
ncvTest(modRM)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.2490362, Df = 1, p = 0.61775spreadLevelPlot(modRM)
##
## Suggested power transformation: 1.818741Modelo con interacción entre Age y Glucose
Para denotar interacción entre variables se usa el símbolo ( : ) Para incluir las variables solas y su interacción se utiliza el símbolo ( * )
modRM <- lm(BMI ~ Cholesterol + Glucose:Age,
data = reg.multiple)
summary(modRM)##
## Call:
## lm(formula = BMI ~ Cholesterol + Glucose:Age, data = reg.multiple)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.5225 -4.6283 -0.8796 3.6945 18.8840
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.879e+01 4.645e+00 4.044 0.000165 ***
## Cholesterol 4.906e-02 2.441e-02 2.010 0.049322 *
## Glucose:Age 3.662e-04 2.057e-04 1.780 0.080625 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.274 on 55 degrees of freedom
## Multiple R-squared: 0.1292, Adjusted R-squared: 0.09751
## F-statistic: 4.079 on 2 and 55 DF, p-value: 0.02229aic <- AIC(modRM)
sprintf("AIC = %.2f", aic)## [1] "AIC = 382.54"# prueba de homocedasticidad
ncvTest(modRM)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.07666205, Df = 1, p = 0.78187spreadLevelPlot(modRM)
##
## Suggested power transformation: 0.9950719Selección automática de modelo - método “stepwise”
Existen métodos para seleccionar automáticamente el mejor modelo, a base de estadísticos indicadores, y que conlleva un procedimiento iterativo. Uno de estos procedimientos es conocido como ‘stepwise’ (por pasos), y aunque no es el más aceptado en la actualidad, ha sido muy usado y es una buena manera de ilustrar el procedimiento, usando nuestros datos.
En este procedimiento el proceso de selección se basa en mantener el modelo con el menor valor del estadístico AIC (Akaike Information Criterion), que indica el modelo con la menor pérdida de información y mayor simplicidad. En el proceso se parte de un modelo nulo (no efecto de predictores) y hasta un modelo muy complejo, incluyendo interacciones. Las variables se incluyen y se quitan, y cada vez se recalcula AIC, hasta obtener el modelo que mantiene el mínimo valor de AIC.
#formulación de un modelo nulo y un modelo completo
modNulo <- lm(BMI ~ 1, data = reg.multiple)
modFull <- lm(BMI ~ Cholesterol*Glucose + Age*Cholesterol + Age*Glucose,
data = reg.multiple)
#procedimiento stepwise
bmistep <- step(modNulo,
scope = list(lower=modNulo, upper=modFull,
direction="both"))## Start: AIC=219.97
## BMI ~ 1
##
## Df Sum of Sq RSS AIC
## + Cholesterol 1 196.462 2289.8 217.19
## + Glucose 1 121.605 2364.6 219.06
## + Age 1 86.322 2399.9 219.92
## <none> 2486.2 219.97
##
## Step: AIC=217.19
## BMI ~ Cholesterol
##
## Df Sum of Sq RSS AIC
## + Glucose 1 98.063 2191.7 216.66
## <none> 2289.8 217.19
## + Age 1 49.203 2240.6 217.93
## - Cholesterol 1 196.462 2486.2 219.97
##
## Step: AIC=216.66
## BMI ~ Cholesterol + Glucose
##
## Df Sum of Sq RSS AIC
## <none> 2191.7 216.66
## + Cholesterol:Glucose 1 55.387 2136.3 217.17
## - Glucose 1 98.063 2289.8 217.19
## + Age 1 21.337 2170.3 218.09
## - Cholesterol 1 172.920 2364.6 219.06summary(bmistep)##
## Call:
## lm(formula = BMI ~ Cholesterol + Glucose, data = reg.multiple)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.125 -4.748 -1.134 3.321 18.832
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.95236 4.80502 3.736 0.000446 ***
## Cholesterol 0.05098 0.02447 2.083 0.041902 *
## Glucose 0.02254 0.01437 1.569 0.122449
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.313 on 55 degrees of freedom
## Multiple R-squared: 0.1185, Adjusted R-squared: 0.08641
## F-statistic: 3.696 on 2 and 55 DF, p-value: 0.0312Comparación de modelos
Gráfica de coeficientes
Una manera de comparar visualmente modelos (en realidad sus coeficientes) es usar el paquete coefplot, en conjunto con ggplot2, para crear una gráfica de los coeficientes estimados de cada variable (sola o de interacción), en cada modelo y detectar los que son diferentes de 0, y los modelos que los contienen.
library(ggplot2)
library(coefplot)##
## Attaching package: 'coefplot'## The following object is masked from 'package:Rmisc':
##
## multiplot#cálculo para todos los modelos
modbmi1 <- lm(BMI ~ Age + Cholesterol + Glucose, data=reg.multiple)
modbmi2 <- lm(BMI ~ Age*Glucose + Cholesterol*Glucose + Age*Cholesterol, data=reg.multiple)
modbmi3 <- lm(BMI ~ Cholesterol + Age:Glucose, data=reg.multiple)
modbmi4 <- lm(BMI ~ Cholesterol + Glucose, data=reg.multiple)
modbmi5 <- lm(BMI ~ Cholesterol, data = reg.multiple)
#comparando coeficientes de todos los modelos
multiplot(modbmi1, modbmi2, modbmi3, modbmi4, modbmi5, pointSize = 2, intercept=FALSE)
Selección de modelo usando R-cuadrado ajustado y Mallow’s Cp para mejores modelos
El estadístico \(R^2\) es la cantidad de varianza en la respuesta (variable dependiente) producido por las variables predictoras, mediante el modelo, por lo tanto constituye una buena manera de medir la capacidad del modelo para “explicar” los datos. Sin embargo, en un modelo de regresión múltiple, el \(R^2\) aumenta al aumentar el número de predictores, lo cual conlleva a sobrestimar la “calidad” del modelo, con un número excesivo de variables. Usando el \(R^2ajustado\) se toma en cuenta el número de parámetros en el modelo, por lo tanto es una medida más realista del ajuste al modelo.
El estadístico de Mallow, \(C_p\), es otro indicador para seleccionar el mejor modelo en una regresión múltiple; funciona de manera similar al AIC, y sirve para evitar incluir parámetros en exceso en el modelo. La regla general es escoger el modelo con el número de parámetros, en el cual el valor de \(C_p\) sea cercano (pero menor) al número de parámetros más 1.
library(leaps)
modSS <- regsubsets(BMI ~ Age*Cholesterol + Cholesterol*Glucose + Age*Glucose, data = reg.multiple, nbest = 3, intercept = TRUE)
# gráfica para R^2 ajustado
plot(modSS, scale="adjr2")
# gráfica para Cp
plot(modSS, scale="Cp")
# otra forma de visualizar Cp
library(car)
## Mallow Cp
mallowCp <-
subsets(modSS, statistic="cp", legend = FALSE, min.size = 1, main = "Mallow Cp")
abline(a = 1, b = 1, lty = 2)
EJERCICIO
Usar los datos deathsmall_cities.xlsx para buscar un modelo de regresión múltiple entre la variable dependiente death1K (muertes anuales por cada 1000 habitantes) y las otras variables del archivo.
=================================
DIA 2
ANOVA
Cuando necesitamos comparar las medias de más de dos muestras, la metodología estadística utilizada comúnmente para estos casos, es el análisis de varianza.
Modelo de Una-Vía (One-Way)
Tenemos el siguiente modelo para describir el valor de cualquier observación, proveniente de varios grupos (\(n_i, i=1...k\)), cada uno con \(n_{ij}\) observaciones:
\[y_{ij} = \mu + \alpha_i + e_{ij}\] donde:
- \(\mu\) representa la media de todas las mediciones de todos los grupos.
- \(\alpha_i\) representa la diferencia entre la media del grupo \(n_i\) y la media de todos los grupos.
- \(e_{ij}\) representa el error aleatorio alrededor de los valores \(\mu + \alpha_i\), de una observación del grupo \(n_i\).
La hipótesis nula para el análisis de varianza de una vía, sería que las medias de los grupos son las mismas (e iguales a la media global, \(\mu\)), para esto todos los valores de \(\alpha_i\) deben ser igual a 0. La hipótesis alterna sería entonces, que al menos un \(\alpha_i\) debe ser diferente de 0.
Cálculo de la variabilidad en forma de Suma de Cuadrados (SS)
Para propósitos de cálculo del estadístico a ser usado para la prueba de hipótesis (F), se calculan los siguientes estimadores de la variabilidad total, dentro de grupos, y entre grupos:
Suma Total de Cuadrados (Total SS)
\[\sum_{i=1}^{k} \sum_{j=1}^{n_i}(y_{ij} - \bar{\bar y})^2\]
Suma de Cuadrados Dentro de Grupos (Within SS)
\[\sum_{i=1}^{k} \sum_{j=1}^{n_i}(y_{ij} - \bar y_i)^2\]
Suma de Cuadrados Entre Grupos (Between SS)
\[\sum_{i=1}^{k} \sum_{j=1}^{n_i}(\bar y_i - \bar{\bar y})^2\]
Ejemplo de tabla de resultados de prueba de hipótesis ANOVA de una vía
Ejemplo 2 y Datos
(tomado de Zar, 2014)
En una investigación, 19 cerdos jóvenes fueron asignados, al azar, a
cuatro grupos experimentales. Cada grupo se alimentó con una dieta
diferente (D1, D2, D3, D4). Luego de ser criados hasta adultos, se midió
la masa corporal (kg) de cada animal. Queremos saber si la masa corporal
resultó ser igual (\(H_0\)) para las
cuatro dietas.
\[H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4\] Vamos a crear una tabla de datos para el ANOVA a partir de vectores; esto es práctico si hay pocos datos y no se quiere depender de un archivo externo:
# creación de vectores de dietas
D1 <- c(60.8,67,65,68.6,61.7)
D2 <- c(68.7,67.7,75,73.3,71.8)
D3 <- c(69.6,77.1,75.2,71.5)
D4 <- c(61.9,64.2,63.1,66.7,60.3)
#data frame con factor Dieta
df.cerdos <- data.frame(Dieta=c(rep("D1", times=length(D1)),
rep("D2", times=length(D2)),
rep("D3", times=length(D3)),
rep("D4", times=length(D4))),
masa=c(D1, D2, D3, D4))
df.cerdos## Dieta masa
## 1 D1 60.8
## 2 D1 67.0
## 3 D1 65.0
## 4 D1 68.6
## 5 D1 61.7
## 6 D2 68.7
## 7 D2 67.7
## 8 D2 75.0
## 9 D2 73.3
## 10 D2 71.8
## 11 D3 69.6
## 12 D3 77.1
## 13 D3 75.2
## 14 D3 71.5
## 15 D4 61.9
## 16 D4 64.2
## 17 D4 63.1
## 18 D4 66.7
## 19 D4 60.3Visualización de los datos
# creating box-plot graph for each group
library(ggplot2)
ggplot(df.cerdos, aes(x=Dieta, y=masa)) + geom_boxplot(fill="cornflowerblue") +
labs(x = "Dietas", y = "Masa corporal, kg")
Función aov para ANOVA
La función básica aov de R calcula los estadísticos del análisis de varianza:
analisis <- aov(masa ~ Dieta, df.cerdos)
summary(analisis)## Df Sum Sq Mean Sq F value Pr(>F)
## Dieta 3 338.9 112.98 12.04 0.000283 ***
## Residuals 15 140.8 9.38
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La prueba de hipótesis del ANOVA, usando el estadístico F, nos indica que la probabilidad del error tipo I es bastante baja (<< 0.05), por lo cual podemos rechazar la hipótesis nula y decir que al menos una dieta produce una mayor ganancia en masa de los cerdos.
Supuestos del ANOVA
La prueba de ANOVA, al igual que la prueba t, asume una distribución normal de los datos, y homocedasticidad (homogeneidad de varianzas). Por lo tanto, antes de continuar con otras pruebas es importante primero determinar si los datos cumplen con estos supuestos.
Prueba de normalidad
Examinaremos gráficamente si los efectos (en realidad los residuales) tienen una distribución normal, usando una gráfica de cuantiles-cuantiles (‘Q-Q plot’):
library(EnvStats)
# response <- lm(masa ~ Dieta, data = df.cerdos)
EnvStats::qqPlot(residuals(analisis), add.line = TRUE, points.col = "blue", line.col = "red")
A continuación la prueba de Shapiro-Wilk para
normalidad, con una hipótesis nula de que los valores se ajustan a la
distribución normal:
shapiro.test(residuals(analisis))##
## Shapiro-Wilk normality test
##
## data: residuals(analisis)
## W = 0.91042, p-value = 0.07534No podemos rechazar la hipótesis nula con esta probabilidad (error tipo I), y aceptamos que los datos tienen una distribución normal.
Para la homogeneidad de varianza realizamos la prueba de Bartlett, que tiene una hipótesis nula de varianzas iguales:
bartlett.test(masa ~ Dieta, data=df.cerdos)##
## Bartlett test of homogeneity of variances
##
## data: masa by Dieta
## Bartlett's K-squared = 0.47515, df = 3, p-value = 0.9243No podemos rechazar la hipótesis nula con esta probabilidad (error tipo I), y aceptamos que los datos por grupo tienen homogeneidad de varianza.
Pruebas post-hoc
El análisis de varianza nos indica, para cierto nivel de probabilidad, si un factor (o más de uno) produce desviaciones de un valor medio global; sin embargo no nos indica cuál o cuáles niveles del factor son los determinantes en la desviación.
Gráfica con intervalos de confianza
Mediante la siguiente gráfica, de medias e intervalos de confianza (95%), podemos hipotetizar cuáles dietas son diferentes entre sí y cuáles no.
library(plyr)
library(ggplot2)
ic.dietas <- ddply(df.cerdos, "Dieta", plyr::summarize,
masa.mean=mean(masa), masa.sd=sd(masa),
Lenght=NROW(masa),
tfrac=qt(p=.95, df=Lenght-1),
Lower=masa.mean - tfrac*masa.sd/sqrt(Lenght),
Upper=masa.mean + tfrac*masa.sd/sqrt(Lenght)
)
ggplot(ic.dietas, aes(x=masa.mean, y=Dieta)) + geom_point() + geom_errorbarh(aes(xmin=Lower, xmax=Upper), height=.3)
Para poder determinar si las diferencias específicas entre los niveles de los factores son significativas, debemos utilizar las llamadas pruebas post-hoc o procedimientos de comparaciones múltiples (MCP).
Utilizaremos la prueba de Tukey para probar cuáles son las dietas del Ejemplo 2 que se diferencian entre si.
library(agricolae)
posthoc <- TukeyHSD(analisis)
posthoc## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = masa ~ Dieta, data = df.cerdos)
##
## $Dieta
## diff lwr upr p adj
## D2-D1 6.68 1.096263 12.263737 0.0168421
## D3-D1 8.73 2.807553 14.652447 0.0034914
## D4-D1 -1.38 -6.963737 4.203737 0.8906642
## D3-D2 2.05 -3.872447 7.972447 0.7530266
## D4-D2 -8.06 -13.643737 -2.476263 0.0041505
## D4-D3 -10.11 -16.032447 -4.187553 0.0009497par(cex = 0.6)
plot(posthoc)
ANOVA con repeticiones (Within-Subject)
En el ANOVA con repeticiones, los sujetos son medidos más de una vez.
Ejemplo 3 y Datos
Los datos de esta investigación provienen de un experimento con ratones a los que se le injertó un tejido con carcinoma humano. Tres grupos de ratones se sometieron a tres tratamientos: control, y dos drogas (Docetaxel y Paclitaxel), y midieron el tamaño de los tumores (\(mm^3\)) seis veces en cada ratón (0, 2, 4, 7, 9, y 11 días). Queremos probar: (a) si hubo un efecto reductor de los tumores por efecto de las drogas (‘Between-Groups’), (b) si el efecto depende del tiempo (‘Within-Subjects’).
Los datos deben tener un arreglo de sujetos x factores-variables.
# datos desde Excel - tomando la hoja 'tumor'
library(readxl)
tumor <- read_excel("Data/tumor_drug_repeat.xlsx",
sheet = "tumor")
head(tumor, 10)## # A tibble: 10 × 4
## mouse drugT time tsize
## <chr> <chr> <chr> <dbl>
## 1 1_1 control day00 296.
## 2 1_1 control day02 343.
## 3 1_1 control day04 354.
## 4 1_1 control day07 396.
## 5 1_1 control day09 387.
## 6 1_1 control day11 432.
## 7 1_2 control day00 190.
## 8 1_2 control day02 226.
## 9 1_2 control day04 212.
## 10 1_2 control day07 284.tail(tumor)## # A tibble: 6 × 4
## mouse drugT time tsize
## <chr> <chr> <chr> <dbl>
## 1 3_12 Paclitaxel day00 221.
## 2 3_12 Paclitaxel day02 334.
## 3 3_12 Paclitaxel day04 231.
## 4 3_12 Paclitaxel day07 177.
## 5 3_12 Paclitaxel day09 172.
## 6 3_12 Paclitaxel day11 326.Gráfica box-plot para visualizar los datos
#gráfica de grupos
boxplot(tsize ~ drugT*time, data = tumor, col = (c("green", "blue", "red")),
las=2, cex.axis=0.70,
main = "Tamaño de tumor x tratamiento",
ylab = "Volumen del tumor, mm^3",
cex.lab=0.75)
# contribución de Robert
library(ggplot2)
ggplot(tumor,aes(x=time,y=tsize,))+
geom_boxplot(aes(fill=factor(drugT)))+
facet_wrap(~drugT)+
theme_minimal()+
ggtitle("tumor size through time and treatment") +
theme(axis.text.x=element_text(angle = 90, vjust = 0.5))
aov para ANOVA con repeticiones
#anova con repeticiones
anovarep <- aov(tsize ~ drugT*time + Error(mouse/(time)), tumor)
summary(anovarep)##
## Error: mouse
## Df Sum Sq Mean Sq F value Pr(>F)
## drugT 2 64468 32234 2.325 0.116
## Residuals 29 402090 13865
##
## Error: mouse:time
## Df Sum Sq Mean Sq F value Pr(>F)
## time 5 21095 4219 1.284 0.27403
## drugT:time 10 80822 8082 2.459 0.00959 **
## Residuals 145 476515 3286
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pruebas de supuestos
# prueba de normalidad de residuales
shapiro.test(residuals(anovarep[[2]]))##
## Shapiro-Wilk normality test
##
## data: residuals(anovarep[[2]])
## W = 0.95444, p-value = 0.2067# prueba de homogeneidad de varianza entre grupos
bartlett.test(tsize ~ drugT, data=tumor)##
## Bartlett test of homogeneity of variances
##
## data: tsize by drugT
## Bartlett's K-squared = 7.9399, df = 2, p-value = 0.01887bartlett.test(tsize ~ time, data=tumor)##
## Bartlett test of homogeneity of variances
##
## data: tsize by time
## Bartlett's K-squared = 36.351, df = 5, p-value = 8.078e-07ANOVA con dos factores (‘two-way ANOVA’)
En el modelo anterior de ANOVA consideramos el efecto de un solo factor (dieta), vamos a estudiar el caso cuando queremos probar si dos factores son determinantes en las diferencias entre las medias de dos o más grupos (también llamados celdas).
Las hipótesis en este caso son más de una:
- hipótesis acerca del efecto del factor A
- hipótesis acerca del efecto del factor B
- hipótesis acerca del efecto de la interacción entre los factores
Ejemplo 4 y Datos
(tomado de Kabacoff, 2015)
Estaremos utilizando datos contenidos en el sistema base de R: ToothGrowth
TG <- ToothGrowth
head(TG)## len supp dose
## 1 4.2 VC 0.5
## 2 11.5 VC 0.5
## 3 7.3 VC 0.5
## 4 5.8 VC 0.5
## 5 6.4 VC 0.5
## 6 10.0 VC 0.5tail(TG)## len supp dose
## 55 24.8 OJ 2
## 56 30.9 OJ 2
## 57 26.4 OJ 2
## 58 27.3 OJ 2
## 59 29.4 OJ 2
## 60 23.0 OJ 2Conejillos de India (n = 60) fueron asignados al azar (10 por grupo, diseño balanceado) a tratamientos combinados (grupos) de tres dosis de ácido ascórbico (dose: 0.5, 1, o 2 mg) y dos tipos de bebidas (supp: jugo de naranja, OJ o bebida con vitamina C, VC). La variable respuesta, luego de un tiempo con los tratamientos, fue la longitud de los dientes frontales (len: mm). El modelo a usar es: \[len = interacción(supp\ x\ dosis)\]
library(gplots) #para la gráfica
#tabla tratamientos
table(TG$supp,TG$dose)##
## 0.5 1 2
## OJ 10 10 10
## VC 10 10 10#medias de los grupos
medias.grupos <- aggregate(TG$len, by=list(TG$supp,TG$dose), FUN=mean)
medias.grupos <- setNames(medias.grupos, c("Bebida", "Dosis, mg", "Media Long. dientes, mm"))
medias.grupos## Bebida Dosis, mg Media Long. dientes, mm
## 1 OJ 0.5 13.23
## 2 VC 0.5 7.98
## 3 OJ 1.0 22.70
## 4 VC 1.0 16.77
## 5 OJ 2.0 26.06
## 6 VC 2.0 26.14#desviación estándar
sd.grupos <- aggregate(TG$len, by=list(TG$supp,TG$dose), FUN=sd)
sd.grupos <- setNames(sd.grupos, c("Bebida", "Dosis, mg", "DE Long. dientes, mm"))
sd.grupos## Bebida Dosis, mg DE Long. dientes, mm
## 1 OJ 0.5 4.459709
## 2 VC 0.5 2.746634
## 3 OJ 1.0 3.910953
## 4 VC 1.0 2.515309
## 5 OJ 2.0 2.655058
## 6 VC 2.0 4.797731#convertir dosis (numérico) a factor
TG$dose <- factor(TG$dose)
TG## len supp dose
## 1 4.2 VC 0.5
## 2 11.5 VC 0.5
## 3 7.3 VC 0.5
## 4 5.8 VC 0.5
## 5 6.4 VC 0.5
## 6 10.0 VC 0.5
## 7 11.2 VC 0.5
## 8 11.2 VC 0.5
## 9 5.2 VC 0.5
## 10 7.0 VC 0.5
## 11 16.5 VC 1
## 12 16.5 VC 1
## 13 15.2 VC 1
## 14 17.3 VC 1
## 15 22.5 VC 1
## 16 17.3 VC 1
## 17 13.6 VC 1
## 18 14.5 VC 1
## 19 18.8 VC 1
## 20 15.5 VC 1
## 21 23.6 VC 2
## 22 18.5 VC 2
## 23 33.9 VC 2
## 24 25.5 VC 2
## 25 26.4 VC 2
## 26 32.5 VC 2
## 27 26.7 VC 2
## 28 21.5 VC 2
## 29 23.3 VC 2
## 30 29.5 VC 2
## 31 15.2 OJ 0.5
## 32 21.5 OJ 0.5
## 33 17.6 OJ 0.5
## 34 9.7 OJ 0.5
## 35 14.5 OJ 0.5
## 36 10.0 OJ 0.5
## 37 8.2 OJ 0.5
## 38 9.4 OJ 0.5
## 39 16.5 OJ 0.5
## 40 9.7 OJ 0.5
## 41 19.7 OJ 1
## 42 23.3 OJ 1
## 43 23.6 OJ 1
## 44 26.4 OJ 1
## 45 20.0 OJ 1
## 46 25.2 OJ 1
## 47 25.8 OJ 1
## 48 21.2 OJ 1
## 49 14.5 OJ 1
## 50 27.3 OJ 1
## 51 25.5 OJ 2
## 52 26.4 OJ 2
## 53 22.4 OJ 2
## 54 24.5 OJ 2
## 55 24.8 OJ 2
## 56 30.9 OJ 2
## 57 26.4 OJ 2
## 58 27.3 OJ 2
## 59 29.4 OJ 2
## 60 23.0 OJ 2#anova dos factores
anova2 <- aov(TG$len ~ TG$supp*TG$dose)
summary(anova2)## Df Sum Sq Mean Sq F value Pr(>F)
## TG$supp 1 205.4 205.4 15.572 0.000231 ***
## TG$dose 2 2426.4 1213.2 92.000 < 2e-16 ***
## TG$supp:TG$dose 2 108.3 54.2 4.107 0.021860 *
## Residuals 54 712.1 13.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#gráfica de medias y intervalo confianza 95%
plotmeans(TG$len ~ interaction(TG$supp,TG$dose,sep = " "),
connect = list(c(1,3,5),c(2,4,6)),
col = c("orange", "darkgreen"),
main = "Gráfica de Medias de Grupos con IC 95%",
xlab = "Combinación de dosis y bebida",
ylab = "Longitud de dientes frontales, mm")
ANCOVA
Aunque en una prueba ANOVA de una vía estamos asumiendo que el único o mayor efecto sobre la variable respuesta, es el factor considerado o tratamiento, otras variables pueden actuar como pseudofactores, produciendo un efecto conjunto con el factor considerado, y las llamamos covariables. Vamos a examinar este tipo de análisis para un ANOVA de una vía (un factor), con una covariable. A este tipo de análisis se le conoce como análisis de covarianza (ANCOVA).
Ejemplo 5 y Datos
(tomado de Kabacoff, 2015)
Ratonas preñadas fueron divididas en cuatro grupos de tratamientos, con diferentes dosis de una droga (dose: 0, 5, 50, o 500). El peso promedio de la camada (weigth) era la variable dependiente; el tiempo de gestación (gesttime) fue la covariable.
#paquete requerido con los datos
library(multcomp)
#datos litter y tabla de n por tratamientos
data(litter)
head(litter)## dose weight gesttime number
## 1 0 28.05 22.5 15
## 2 0 33.33 22.5 14
## 3 0 36.37 22.0 14
## 4 0 35.52 22.0 13
## 5 0 36.77 21.5 15
## 6 0 29.60 23.0 5table(litter$dose)##
## 0 5 50 500
## 20 19 18 17#media por grupo de tratamiento
medias.trat <- aggregate(litter$weight, by=list(litter$dose), FUN=mean)
medias.trat <- setNames(medias.trat, c("Dosis", "Masa, g"))
medias.trat## Dosis Masa, g
## 1 0 32.30850
## 2 5 29.30842
## 3 50 29.86611
## 4 500 29.64647#ANCOVA
ancova <- aov(weight ~ gesttime + dose, data=litter)
summary(ancova)## Df Sum Sq Mean Sq F value Pr(>F)
## gesttime 1 134.3 134.30 8.049 0.00597 **
## dose 3 137.1 45.71 2.739 0.04988 *
## Residuals 69 1151.3 16.69
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Los resultados muestran un efecto significativo del tiempo de gestación, y marginalmente significativo de la dosis. Esto quiere decir que las medias de masa calculadas para los grupos tratados, están afectadas por el tiempo de gestación; es posible remover este efecto, usando la función effect:
library(effects)
effect("dose", ancova)##
## dose effect
## dose
## 0 5 50 500
## 32.35367 28.87672 30.56614 29.33460MANOVA - Análisis de Varianza Multivariado
Si tenemos más de una variable respuesta, se puede probar el efecto de factores de manera simultánea en todas las variables, usando un MANOVA.
Ejemplo 6 y Datos
En este ejemplo se investiga el efecto de la posición de caja de cereales, de diversas marcas, a diferentes niveles de almacenamiento, en el contenido de calorías, grasa y azúcares, que vendrían siendo las variables dependientes; este sería una MANOVA de un factor, con tres niveles: 1 (parte baja), 2 (parte media), y 3 (parte alta).
library(MASS)
# datos en R básico
data(UScereal) # cargar los datos
head(UScereal)## mfr calories protein fat sodium fibre
## 100% Bran N 212.1212 12.121212 3.030303 393.9394 30.303030
## All-Bran K 212.1212 12.121212 3.030303 787.8788 27.272727
## All-Bran with Extra Fiber K 100.0000 8.000000 0.000000 280.0000 28.000000
## Apple Cinnamon Cheerios G 146.6667 2.666667 2.666667 240.0000 2.000000
## Apple Jacks K 110.0000 2.000000 0.000000 125.0000 1.000000
## Basic 4 G 173.3333 4.000000 2.666667 280.0000 2.666667
## carbo sugars shelf potassium vitamins
## 100% Bran 15.15152 18.18182 3 848.48485 enriched
## All-Bran 21.21212 15.15151 3 969.69697 enriched
## All-Bran with Extra Fiber 16.00000 0.00000 3 660.00000 enriched
## Apple Cinnamon Cheerios 14.00000 13.33333 1 93.33333 enriched
## Apple Jacks 11.00000 14.00000 2 30.00000 enriched
## Basic 4 24.00000 10.66667 3 133.33333 enriched# transformar una variable numérica a factor
shelf <- factor(UScereal$shelf)
# selección de variables dependientes en una matriz
y <- cbind(UScereal$calories, UScereal$fat, UScereal$sugars)
head(y)## [,1] [,2] [,3]
## [1,] 212.1212 3.030303 18.18182
## [2,] 212.1212 3.030303 15.15151
## [3,] 100.0000 0.000000 0.00000
## [4,] 146.6667 2.666667 13.33333
## [5,] 110.0000 0.000000 14.00000
## [6,] 173.3333 2.666667 10.66667# media de las variables dependientes por grupo (shelf)
aggregate(y, by=list(shelf), FUN=mean)## Group.1 V1 V2 V3
## 1 1 119.4774 0.6621338 6.295493
## 2 2 129.8162 1.3413488 12.507670
## 3 3 180.1466 1.9449071 10.856821Uso de la función manova
fit <- manova(y ~ shelf)
summary(fit)## Df Pillai approx F num Df den Df Pr(>F)
## shelf 2 0.4021 5.1167 6 122 0.0001015 ***
## Residuals 62
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como obtuvimos un efecto significativo en la prueba, podemos usar los resultados de los ANOVAs de cada variable dependiente:
summary.aov(fit)## Response 1 :
## Df Sum Sq Mean Sq F value Pr(>F)
## shelf 2 50435 25217.6 7.8623 0.0009054 ***
## Residuals 62 198860 3207.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response 2 :
## Df Sum Sq Mean Sq F value Pr(>F)
## shelf 2 18.44 9.2199 3.6828 0.03081 *
## Residuals 62 155.22 2.5035
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response 3 :
## Df Sum Sq Mean Sq F value Pr(>F)
## shelf 2 381.33 190.667 6.5752 0.002572 **
## Residuals 62 1797.87 28.998
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1===============================
DIA 3
GLM
En los días 1 y 2 del taller, examinamos modelos lineales que pueden ser usados para probar o predecir la respuesta de una variable respuesta con distribución normal, a partir de variables predictoras continuas y/o categóricas. Pero hay muchas situaciones en las cuales la variable respuesta o dependiente no está distribuida normalmente, o hasta no es continua. Para estos casos podemos usar lo que se conoce como ‘Generalized Linear Models’ (GLM). Entre estos tenemos:
- Cuando la variable dependiente es categórica de tipo binaria (0, 1; si, no; vivo, muerto) o politoma (malo, regular, bueno), se puede usar la regresión logística.
- Cuando la variable dependiente es un conteo (número de casos por día), los valores se van a limitar a ciertos números enteros, y nunca negaticos; en este caso se utiliza la regresión con distribución Poisson.
Regresión Logística
A continuación voy a mostrar un ejemplo del uso de regresión logística, utilizando la función glm (parte del paquete stats de la instalación básica de R).
Ejemplo 7 y Datos
Los datos que vamos a analizar son de los pasajeros del viaje del Titanic que naufragó. Estos datos muestran los detalles de cada pasajero (variables predictoras) y si cada pasajero sobrevivió o no (variable dependiente binaria). El objetivo es analizar los datos y encontrar un modelo de regresión logística que nos indique la(s) variable(s) que mejor predicen la supervivencia de los pasajeros. Las variables a utilizar son las siguientes:
- Survived: 0 no, 1 si
- Pclass: categoría del camarote (1, primera; 2, segunda; 3, tercera)
- Sex: male o female
- Age: edad
- Nfamily: número de familiares viajando, se va a construir sumando SibSp y Parch.
Cargamos los datos a partir de una hoja de cálculo externa:
library(readxl)
Titanic <- read_excel("Data/Titanic.xlsx")
head(Titanic)## # A tibble: 6 × 12
## PassengerId Survived Pclass Name Sex Age SibSp Parch Ticket Fare Cabin
## <dbl> <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 1 0 3 Braund… male 22 1 0 A/5 2… 7.25 <NA>
## 2 2 1 1 Cuming… fema… 38 1 0 PC 17… 71.3 C85
## 3 3 1 3 Heikki… fema… 26 0 0 STON/… 7.92 <NA>
## 4 4 1 1 Futrel… fema… 35 1 0 113803 53.1 C123
## 5 5 0 3 Allen,… male 35 0 0 373450 8.05 <NA>
## 6 6 0 3 Moran,… male NA 0 0 330877 8.46 <NA>
## # ℹ 1 more variable: Embarked <chr>Creamos la variable Nfamily:
library(dplyr)
Titanic <- Titanic %>% mutate(Nfamily = SibSp + Parch)
head(Titanic)## # A tibble: 6 × 13
## PassengerId Survived Pclass Name Sex Age SibSp Parch Ticket Fare Cabin
## <dbl> <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 1 0 3 Braund… male 22 1 0 A/5 2… 7.25 <NA>
## 2 2 1 1 Cuming… fema… 38 1 0 PC 17… 71.3 C85
## 3 3 1 3 Heikki… fema… 26 0 0 STON/… 7.92 <NA>
## 4 4 1 1 Futrel… fema… 35 1 0 113803 53.1 C123
## 5 5 0 3 Allen,… male 35 0 0 373450 8.05 <NA>
## 6 6 0 3 Moran,… male NA 0 0 330877 8.46 <NA>
## # ℹ 2 more variables: Embarked <chr>, Nfamily <dbl>Gráficas logísticas
Vamos a graficar supervivencia en función de variables continuas:
library(ggplot2)
#gráfica con curva logística
ggplot(Titanic, aes(x=Nfamily, y=Survived, na.rm = TRUE)) +
geom_point() +
geom_smooth(method = "glm",
method.args = list(family = "binomial"),
se = FALSE)## `geom_smooth()` using formula = 'y ~ x'
ggplot(Titanic, aes(x=Age, y=Survived, na.rm = TRUE)) +
geom_point() +
geom_smooth(method = "glm",
method.args = list(family = "binomial"),
se = FALSE)## `geom_smooth()` using formula = 'y ~ x'
Tablas de contraste
Podemos obtener los valores de número de casos en tablas de contraste entre la supervivencia y las variables categóricas:
table(Titanic$Survived, Titanic$Sex)##
## female male
## 0 92 734
## 1 374 109table(Titanic$Survived, Titanic$Pclass)##
## 1 2 3
## 0 137 160 529
## 1 186 117 180Uso de la función glm para probar modelos de regresión logística
#modelo logístico:
reglog <- glm(Survived ~ Age + Sex + Pclass + Nfamily, data = Titanic,
family=binomial(logit))
#resumen de los resultados:
summary(reglog) ##
## Call:
## glm(formula = Survived ~ Age + Sex + Pclass + Nfamily, family = binomial(logit),
## data = Titanic)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.418827 0.448752 9.847 < 2e-16 ***
## Age -0.013811 0.006591 -2.095 0.03615 *
## Sexmale -3.586491 0.202961 -17.671 < 2e-16 ***
## Pclass -0.997279 0.121686 -8.195 2.5e-16 ***
## Nfamily -0.178687 0.065789 -2.716 0.00661 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1401.72 on 1045 degrees of freedom
## Residual deviance: 816.53 on 1041 degrees of freedom
## (263 observations deleted due to missingness)
## AIC: 826.53
##
## Number of Fisher Scoring iterations: 5Para interpretar mejor las diferencias entre el efecto de las variables predictoras, es mejor usar ‘razones de probabilidades’ (en inglés odd ratios (ORs)), o proporción de casos favorables a desfavorables:
\[OR_i = \frac{\pi}{1-\pi} \qquad \] donde \(\pi\) es la probabilidad de que el evento dependiente (\(Y\)) ocurra.
#obtener los coeficientes como ORs:
exp(cbind(OR = coef(reglog), confint(reglog)))## Waiting for profiling to be done...## OR 2.5 % 97.5 %
## (Intercept) 82.99889307 35.13335849 204.38944482
## Age 0.98628432 0.97349704 0.99901923
## Sexmale 0.02769535 0.01838577 0.04078259
## Pclass 0.36888197 0.28929248 0.46643303
## Nfamily 0.83636725 0.73348487 0.94976234Aprendizaje Automático - Machine Learning
Los procedimientos de regresión producen modelos que se ubican en el campo de lo que se conoce como Aprendizaje Automático (‘Machine Learning’).
A continuación utilizaremos el modelo de regresión logística para ilustrar una forma rudimentaria de proceso de aprendizaje automático.
División de los datos
En primer lugar vamos a separar los datos en dos grupos: uno de aprendizaje (training) y otro de prueba o validación:
library(caret)
# Split the data into training and test set
set.seed(123)
training.samples <- Titanic$Survived %>%
createDataPartition(p = 0.8, list = FALSE)
train.data <- Titanic[training.samples, ]
test.data <- Titanic[-training.samples, ]Modelo usando glm
Obtenemos el modelo utilizando los datos de aprendizaje:
library(MASS)
# Fit the model
model1 <- glm(Survived ~ Age + Sex + Pclass + Nfamily, data = train.data, family = binomial(logit))
# Summarize the final selected model
summary(model1)##
## Call:
## glm(formula = Survived ~ Age + Sex + Pclass + Nfamily, family = binomial(logit),
## data = train.data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.159842 0.488266 8.520 < 2e-16 ***
## Age -0.008157 0.007110 -1.147 0.2513
## Sexmale -3.535175 0.226353 -15.618 < 2e-16 ***
## Pclass -0.977438 0.134437 -7.271 3.58e-13 ***
## Nfamily -0.163385 0.072316 -2.259 0.0239 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1115.56 on 829 degrees of freedom
## Residual deviance: 665.38 on 825 degrees of freedom
## (218 observations deleted due to missingness)
## AIC: 675.38
##
## Number of Fisher Scoring iterations: 5Predicciones usando el 20% de los datos
Ahora usamos los datos de prueba para calcular el nivel de predicción del modelo:
# Make predictions
probabilities <- model1 %>% predict(test.data, type = "response")
predicted.classes <- ifelse(probabilities > 0.5, "1", "0")
# Prediction accuracy
observed.classes <- test.data$Survived
mean(predicted.classes == observed.classes, na.rm = TRUE)## [1] 0.8796296Probando con otro modelo
library(MASS)
# Fit the model
model2 <- glm(Survived ~ Sex + Pclass + Nfamily, data = train.data, family = binomial(logit))
# Summarize the final selected model
summary(model2)##
## Call:
## glm(formula = Survived ~ Sex + Pclass + Nfamily, family = binomial(logit),
## data = train.data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.91000 0.33352 11.723 < 2e-16 ***
## Sexmale -3.68952 0.20838 -17.706 < 2e-16 ***
## Pclass -0.92656 0.11267 -8.223 < 2e-16 ***
## Nfamily -0.22882 0.06069 -3.770 0.000163 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1381.4 on 1047 degrees of freedom
## Residual deviance: 814.1 on 1044 degrees of freedom
## AIC: 822.1
##
## Number of Fisher Scoring iterations: 5# Make predictions
probabilities <- model2 %>% predict(test.data, type = "response")
predicted.classes <- ifelse(probabilities > 0.5, "1", "0")
# Prediction accuracy
observed.classes <- test.data$Survived
mean(predicted.classes == observed.classes, na.rm = TRUE)## [1] 0.8850575Referencias
Kabacoff, R., 2015. R in action: data analysis and graphics with R, Second edition. ed. Manning, Shelter Island.
Lander, J. P., 2014. R for everyone. Pearson Education, Inc., Upper Saddle River, NJ, USA.
Apéndice
https://www.statmethods.net/stats/anova.html
- One Way Anova (Completely Randomized Design)
- fit <- aov(y ~ A, data=mydataframe)
- Randomized Block Design (B is the blocking factor)
- fit <- aov(y ~ A + B, data=mydataframe)
- Two Way Factorial Design
- fit <- aov(y ~ A + B + A:B, data=mydataframe)
- fit <- aov(y ~ A*B, data=mydataframe) # same thing
- Analysis of Covariance
- fit <- aov(y ~ A + x, data=mydataframe)
For within subjects designs, the data frame has to be rearranged so that each measurement on a subject is a separate observation.
- One Within Factor
- fit <- aov(y~A+Error(Subject/A),data=mydataframe)
- Two Within Factors W1 W2, Two Between Factors B1 B2
- fit <- aov(y~(W1W2B1B2)+Error(Subject/(W1W2))+(B1*B2), data=mydataframe)