ANOVA
# Set working directory
library(car)
library(emmeans)
library(rstatix)
library(tidyverse)
library(ggpubr)
options(max.print = 1e+05)Introducción ANOVA
Uno de los procesos más usados para comparar diferentes “tratamientos” sobre un grupo experimental es sin duda el Análisis de Varianza (ANOVA). Desarrollado originalmente por R.A. Fisher y ha tenido una gran influencia en la estadística. El análisis de varianza se usa principalmente para dos propósitos: 1) para estimar y desarrollar prueba de hipótesis de varianzas de una población. 2) Estimar y probar medias poblacionales. En esta introducción resolveremos ejemplos para la estimación de medias poblacionales.
ANOVA Modelos de una via (One-Way ANOVA)
Carne proveniente de caza, entre otras aquellas del venado de cola-blanca, y ardilla gris, son usadas como alimento por familias, cazadores y otros individuos por razones culturales, personales, o de salud. En un estudio de Holben et al. estimó el contenido de selenio de carne de venado cola-blanca libre (venison), y ardilla gris (squirrel) obtenidos de una región de baja concentración de selenio en EEUU. Estos contenidos de selenio fueron también comparados con concentraciones de carne de ganado producido en la misma región (RRB) y en una región externa (NRRB). Nos interesa saber si la concentración de selenio \((\mu g/100g)\) es distinta en los grupos de carne. (EXA_C08_S02_01.csv).
La hipótesis nula \(H_0\) de este modelo es que todos los promedios poblacionales (cada promedio de la población de cada caso) son iguales entre ellos \(\mu_1 = \mu_2 = \mu_3 \dots = \mu_n\), y la hipótesis alternativa es que “por lo menos uno de ellos es distinto”.
# Exa_C08_S02_01 One-way
Exa8.2.1 <- read_csv(file = "~/Dropbox/GitHub/ProbEstad/DataSets/ch08_all/EXA_C08_S02_01mod.csv",
show_col_types = FALSE)
head(Exa8.2.1)# A tibble: 6 × 2
Group Selenium
<chr> <dbl>
1 VEN 26.7
2 VEN 28.6
3 VEN 29.7
4 VEN 27.0
5 VEN 11.0
6 VEN 22.0
#
ggplot(data = Exa8.2.1, aes(x = Group, y = Selenium)) + geom_boxplot() + geom_jitter(width = 0.1,
alpha = 0.6, aes(color = Group))
# one way ANOVA estimated as a lineas model
Exa8.2.1_lm <- lm(Selenium ~ Group, data = Exa8.2.1)
(anova_Exa8.2.1_lm <- anova(Exa8.2.1_lm))Analysis of Variance Table
Response: Selenium
Df Sum Sq Mean Sq F value Pr(>F)
Group 3 18935 6311.7 22.614 5.345e-12 ***
Residuals 140 39074 279.1
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Df Sum Sq Mean Sq F value Pr(>F)
Group 3 18935 6312 22.61 5.34e-12 ***
Residuals 140 39074 279
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Tables of means
Grand mean
35.44736
Group
NRB RRB SQU VEN
62.05 28.52 37.42 25.88
rep 19.00 30.00 53.00 42.00Tables of effects
Group
NRB RRB SQU VEN
26.6 -6.925 1.97 -9.572
rep 19.0 30.000 53.00 42.000 Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Selenium ~ Group, data = Exa8.2.1)
$Group
diff lwr upr p adj
RRB-NRB -33.523982 -46.260203 -20.787762 0.0000000
SQU-NRB -24.629335 -36.244663 -13.014007 0.0000010
VEN-NRB -36.170840 -48.180851 -24.160828 0.0000000
SQU-RRB 8.894648 -1.030122 18.819418 0.0961032
VEN-RRB -2.646857 -13.030768 7.737053 0.9108599
VEN-SQU -11.541505 -20.515363 -2.567647 0.0057777
Más ejemplos
Pacientes de enfermedades reumáticas u osteoporosis generalmente sufren de pérdidas críticas de densidad mineral osea (Bone Mineral Density, BMD). Un medicamento usado para recuperar o prevenir una pérdida mayor de BMD, es el Alendronato. Holcomb y Rothenberg examinaron a 96 mujeres tomando alendronato para determinar si había alguna diferencia en el promedio de cambio en densidad ósea entre cinco clasificaciones diagnostica primarias. El Grupo1 era de pacientes diagnosticados con artritis reumatoide (RA). Grupo2 era de pacientes con una variedad de diagnósticos incluyendo, Lupus, granulomatosis de Wegener y poliarteritis, y otros desordenes vasculares (LUPUS). Grupo3 consistió en pacientes diagnosticados con polimialgia reumática o artritis temporal (PMRTA). El Grupo4 estaba integrado por pacientes con artrosis (OA). Y el Grupo 5 de pacientes con diagnóstico de osteoporosis (O) sin otros desordenes reumáticos. ¿Puede determinar que los promedios de los grupos presentan alguna diferencia? Haga una gráfica de cajas de los grupos.
# Exer. 8.2.2 Bone Mineral Density, BMD Patients suffering from rheumatic
# diseases or osteoporosis often suffer critical loss in bone mineral density
# (BMD).
EXR8.2.2 <- read_csv(file = "~/Dropbox/GitHub/ProbEstad/DataSets/ch08_all/EXR_C08_S02_02.csv",
show_col_types = FALSE)
EXR8.2.2# A tibble: 96 × 2
BMD GROUP
<dbl> <dbl>
1 11.1 1
2 24.4 1
3 10.0 1
4 -3.16 1
5 6.84 1
6 3.32 1
7 1.49 1
8 -1.86 1
9 5.39 1
10 3.87 1
# ℹ 86 more rows [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
[39] 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4
[77] 4 4 4 4 4 4 4 4 4
[ reached 'max' / getOption("max.print") -- omitted 11 entries ]EXR8.2.2_M <- EXR8.2.2 %>%
mutate(GROUP = case_when(GROUP == 1 ~ "RA", GROUP == 2 ~ "LUPUS", GROUP == 3 ~
"PMRTA", GROUP == 4 ~ "OA", GROUP == 5 ~ "O"))
EXR8.2.2_M <- EXR8.2.2_M %>%
mutate(GROUP = GROUP %>%
fct_relevel("O", "RA", "LUPUS", "PMRTA", "OA"))
EXR8.2.2_M %>%
ggplot(aes(x = GROUP, y = BMD)) + geom_boxplot() + geom_jitter(width = 0.1, alpha = 0.6,
aes(color = GROUP))
#
Exr2_2.aov <- aov(BMD ~ GROUP, data = EXR8.2.2_M)
# anova() function to print the ANOVA table
anova(Exr2_2.aov)Analysis of Variance Table
Response: BMD
Df Sum Sq Mean Sq F value Pr(>F)
GROUP 4 355.5 88.864 2.2772 0.06697 .
Residuals 91 3551.1 39.024
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Variance Table
Response: BMD
Df Sum Sq Mean Sq F value Pr(>F)
factor(GROUP) 4 355.5 88.864 2.2772 0.06697 .
Residuals 91 3551.1 39.024
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Call: lm(formula = BMD ~ factor(GROUP), data = EXR8.2.2_M)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.672 1.975 4.896 4.22e-06 ***
factor(GROUP)RA -5.202 2.226 -2.336 0.02166 *
factor(GROUP)LUPUS -5.094 2.870 -1.775 0.07929 .
factor(GROUP)PMRTA -7.491 2.518 -2.975 0.00375 **
factor(GROUP)OA -4.467 2.351 -1.900 0.06061 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard deviation: 6.247 on 91 degrees of freedom
Multiple R-squared: 0.09099
F-statistic: 2.277 on 4 and 91 DF, p-value: 0.06697
AIC BIC
631.06 646.45 - Ilich-Ernst y colaboradores estudiaron el consumo de calcio en la dieta de un grupo de 113 mujeres adultas sanas de edades entre 20 y 88 años. Los investigadores segregaron a los sujetos de estudio por grupos de edad de la siguiente manera: Grupo A; 20.0 – 45.9 años; grupo B; 46.00 – 55. 9 años; grupo C; 56.0 – 65.9 años; y grupo D; de más de 66 años. El consumo de calcio estuvo medido en mg/día. Los datos están listados en el archivo EXR_C08_S02_03.csv. Siguiendo un ANOVA, ¿se puede concluir que hay una diferencia en los promedios de las poblaciones? Haga un HSD Tukey para estimar las diferencias entre las distintas poblaciones, sea Alpha = .05. Grafique los datos usando una gráfica de cajas y explique los resultados.
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # Exr
# 8.2.3
Exr8.2.3 <- read_csv(file = "~/Dropbox/GitHub/ProbEstad/DataSets/ch08_all/EXR_C08_S02_03.csv",
show_col_types = FALSE)
head(Exr8.2.3)# A tibble: 6 × 2
calcium Group
<dbl> <chr>
1 1820 A
2 2588 A
3 2670 A
4 1022 A
5 1555 A
6 222 A 
ggplot(data = Exr8.2.3, aes(x = Group, y = calcium)) + geom_boxplot() + geom_jitter(width = 0.1,
alpha = 0.5, aes(color = Group))
Analysis of Variance Table
Response: calcium
Df Sum Sq Mean Sq F value Pr(>F)
Group 3 5931208 1977069 9.3588 1.476e-05 ***
Residuals 109 23026500 211252
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = calcium ~ Group, data = Exr8.2.3)
$Group
diff lwr upr p adj
B-A -455.72078 -865.7066 -45.73492 0.0230454
C-A -574.53605 -913.5892 -235.48294 0.0001356
D-A -596.63447 -905.3867 -287.88228 0.0000109
C-B -118.81527 -509.0844 271.45386 0.8568796
D-B -140.91369 -505.1676 223.34020 0.7443051
D-C -22.09842 -304.1436 259.94679 0.9969598emmeans del inglés “Estimated Marginal Means” proporciona métodos para obtener medias marginales estimadas (EMMs, también conocidas como medias de mínimos cuadrados) En R la librería emmeans perminte calcular contrastes o combinaciones lineales de dichas medias marginales con varios ajustes por multiplicidad. También se pueden estimar y contrastar pendientes de líneas de tendencia.
$emmeans
Group emmean SE df lower.CL upper.CL
A 1448 98.0 109 1254 1643
B 993 123.0 109 749 1236
C 874 85.3 109 705 1043
D 852 66.3 109 720 983
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
A - B 455.7 157 109 2.900 0.0230
A - C 574.5 130 109 4.421 0.0001
A - D 596.6 118 109 5.042 <0.0001
B - C 118.8 150 109 0.794 0.8569
B - D 140.9 140 109 1.009 0.7443
C - D 22.1 108 109 0.204 0.9970
P value adjustment: tukey method for comparing a family of 4 estimates $emmeans
Group emmean SE df lower.CL upper.CL
A 1448 98.0 109 1254 1643
B 993 123.0 109 749 1236
C 874 85.3 109 705 1043
D 852 66.3 109 720 983
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
B - A -456 157 109 -2.900 0.0128
C - A -575 130 109 -4.421 <0.0001
D - A -597 118 109 -5.042 <0.0001
P value adjustment: dunnettx method for 3 tests Tables of means
Grand mean
991.0177
Group
A B C D
1448 992.6 873.8 851.7
rep 22 14.0 29.0 48.0Tables of effects
Group
A B C D
457.3 1.625 -117.2 -139.3
rep 22.0 14.000 29.0 48.0ANOVA de dos vías
Para el caso de un ANOVA de dos vías, tenemos que la función repuesta responde a dos factores que a su vez tiene dos o más niveles por factor: \(Y \sim factor_1 + factor_2\). En este caso ya se puede presentar un término de interacción entre factores.
Del ejemplo 8.3.1 de Daniel para un ANOVA de dos vías: Un fisioterapeuta quería comparar tres métodos para enseñar a sus pacientes a usar una prótesis. Pensó que la proporción de aprendizage sería diferente dependiendo de la edad de los parcientes y desarrolló un experimento que le permitió tomar en cuenta la edad. Se construye la tabla de los datos del problema.
# Build a table, index by index since Two-Way ANOVA EXAMPLE 8.3.1 has NO-data
# file, therefore we make the example table.
Learn <- tibble(Age = factor(rep(c(1, 2, 3, 4, 5), 3)), Method = factor(rep(1:3,
c(5, 5, 5))), Rate = c(7, 8, 9, 10, 11, 9, 9, 9, 9, 12, 10, 10, 12, 12, 14))
Learn# A tibble: 15 × 3
Age Method Rate
<fct> <fct> <dbl>
1 1 1 7
2 2 1 8
3 3 1 9
4 4 1 10
5 5 1 11
6 1 2 9
7 2 2 9
8 3 2 9
9 4 2 9
10 5 2 12
11 1 3 10
12 2 3 10
13 3 3 12
14 4 3 12
15 5 3 14
Para las gráficas de cajas ejecutadas por un ggplot (o boxplot) es necesario mandar por separado cada grupo en la variable predictora (o independiente), o separar dos páneles, uno por cada grupo.
# For the Age group
ggplot(data = Learn, aes(x = Age, y = Rate)) + geom_boxplot() + geom_jitter(width = 0.1,
alpha = 0.6, aes(color = Method)) + theme_bw()
# For the Rate group
ggplot(data = Learn, aes(x = Method, y = Rate)) + geom_boxplot() + geom_jitter(width = 0.1,
alpha = 0.8, aes(color = Age)) + theme_bw()
Analysis of Variance Table
Response: Rate
Df Sum Sq Mean Sq F value Pr(>F)
Method 2 18.5333 9.2667 21.385 0.0006165 ***
Age 4 24.9333 6.2333 14.385 0.0010017 **
Residuals 8 3.4667 0.4333
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Rate ~ Method + Age, data = Learn)
$Method
diff lwr upr p adj
2-1 0.6 -0.5896489 1.789649 0.3666717
3-1 2.6 1.4103511 3.789649 0.0006358
3-2 2.0 0.8103511 3.189649 0.0034083
$Age
diff lwr upr p adj
2-1 0.3333333 -1.5235390 2.190206 0.9676094
3-1 1.3333333 -0.5235390 3.190206 0.1877558
4-1 1.6666667 -0.1902056 3.523539 0.0810838
5-1 3.6666667 1.8097944 5.523539 0.0009146
3-2 1.0000000 -0.8568723 2.856872 0.4057524
4-2 1.3333333 -0.5235390 3.190206 0.1877558
5-2 3.3333333 1.4764610 5.190206 0.0017351
4-3 0.3333333 -1.5235390 2.190206 0.9676094
5-3 2.3333333 0.4764610 4.190206 0.0154324
5-4 2.0000000 0.1431277 3.856872 0.0348816
Más Ejemplos de ANOVA de dos vías (two-way ANOVA)
En un estudio de efectos pulmonares en conejillos de Indias, Lacroix et al. expuso a 18 conejillos de Indias sensibilizados y 18 no-sensibilizados con albúmina, a tres tratamientos distintos: aire regular, benzaldehido y acetildehido. Al final de la exposición los conejillos de Indias fueron anesteciados y la respuesta a alergias medida por lavado broncoalveolar (BAL). La tabla de resultados se muestra en REV_C08_16 y contiene la cuenta de células alveolares \((\times 10^6)\) por grupo de tratamiento (aire, benzaldehido y acetildehido) y sensibilizados o no con albúmina. Pruebe por las diferencias en a) entre sensibilizados y no sensibilizados por albúmina, b) entre los tres tratamientos. Sea \(\alpha=.05\) en todas las pruebas.
Rev8.16 <- read_csv("~/Dropbox/GitHub/ProbEstad/DataSets/ch08_all/REV_C08_16.csv",
show_col_types = FALSE)
# how is the data ordered?
Rev8.16# A tibble: 36 × 3
Sens Treat Count
<chr> <chr> <dbl>
1 No Act 49.9
2 No Act 50.6
3 No Act 50.4
4 No Act 44.1
5 No Act 36.3
6 No Act 39.2
7 No Air 24.2
8 No Air 24.6
9 No Air 22.6
10 No Air 25.1
# ℹ 26 more rowsRev8.16_fac <- Rev8.16 %>%
mutate(Sens = factor(Sens), Treat = factor(Treat))
Rev8.16_fac <- Rev8.16_fac %>%
mutate(Treat = Treat %>%
fct_relevel("Air", "Act", "Benz"))
boxplot(Count ~ Sens + Treat, data = Rev8.16_fac)
ggplot(data = Rev8.16_fac, aes(x = Treat, y = Count)) + geom_boxplot() + geom_jitter(width = 0.1,
alpha = 0.6, aes(color = Sens))
ggplot(data = Rev8.16_fac, aes(x = Treat, y = Count, fill = Sens)) + geom_boxplot(position = position_dodge(width = 0.6),
alpha = 0.5) + geom_jitter(position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.6), aes(color = Sens), alpha = 0.8, size = 1.4) + labs(x = "Treatment",
y = "Response in Counts", fill = "Sens") + theme_bw()
Analysis of Variance Table
Response: Count
Df Sum Sq Mean Sq F value Pr(>F)
Sens 1 7906.2 7906.2 20.718 7.276e-05 ***
Treat 2 7688.6 3844.3 10.074 0.0004042 ***
Residuals 32 12211.5 381.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova Table (Type II tests)
Response: Count
Sum Sq Df F value Pr(>F)
Sens 7906.2 1 20.718 7.276e-05 ***
Treat 7688.6 2 10.074 0.0004042 ***
Residuals 12211.5 32
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Variance Table
Response: Count
Df Sum Sq Mean Sq F value Pr(>F)
Sens 1 7906.2 7906.2 20.718 7.276e-05 ***
Treat 2 7688.6 3844.3 10.074 0.0004042 ***
Residuals 32 12211.5 381.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1$emmeans
Sens Treat emmean SE df lower.CL upper.CL
No Air 31.4 6.51 32 18.17 44.7
Yes Air 61.1 6.51 32 47.81 74.3
No Act 48.0 6.51 32 34.69 61.2
Yes Act 77.6 6.51 32 64.33 90.9
No Benz 12.2 6.51 32 -1.07 25.5
Yes Benz 41.8 6.51 32 28.56 55.1
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
No Air - Yes Air -29.64 6.51 32 -4.552 0.0009
No Air - No Act -16.52 7.98 32 -2.071 0.3273
No Air - Yes Act -46.16 10.30 32 -4.483 0.0012
No Air - No Benz 19.25 7.98 32 2.413 0.1820
No Air - Yes Benz -10.39 10.30 32 -1.009 0.9114
Yes Air - No Act 13.12 10.30 32 1.275 0.7963
Yes Air - Yes Act -16.52 7.98 32 -2.071 0.3273
Yes Air - No Benz 48.88 10.30 32 4.748 0.0005
Yes Air - Yes Benz 19.25 7.98 32 2.413 0.1820
No Act - Yes Act -29.64 6.51 32 -4.552 0.0009
No Act - No Benz 35.76 7.98 32 4.484 0.0011
No Act - Yes Benz 6.12 10.30 32 0.595 0.9907
Yes Act - No Benz 65.40 10.30 32 6.352 <0.0001
Yes Act - Yes Benz 35.76 7.98 32 4.484 0.0011
[ reached 'max' / getOption("max.print") -- omitted 1 row ]
P value adjustment: tukey method for comparing a family of 6 estimates $emmeans
Sens Treat emmean SE df lower.CL upper.CL
No Air 31.4 6.51 32 18.17 44.7
Yes Air 61.1 6.51 32 47.81 74.3
No Act 48.0 6.51 32 34.69 61.2
Yes Act 77.6 6.51 32 64.33 90.9
No Benz 12.2 6.51 32 -1.07 25.5
Yes Benz 41.8 6.51 32 28.56 55.1
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Yes Air - No Air 29.6 6.51 32 4.552 0.0003
No Act - No Air 16.5 7.98 32 2.071 0.1723
Yes Act - No Air 46.2 10.30 32 4.483 0.0004
No Benz - No Air -19.2 7.98 32 -2.413 0.0868
Yes Benz - No Air 10.4 10.30 32 1.009 0.7331
P value adjustment: dunnettx method for 5 tests # With interactions
Rev8.16_lm_int <- aov(Count ~ Sens * Treat, data = Rev8.16_fac)
# model with interaction
anova(Rev8.16_lm_int)Analysis of Variance Table
Response: Count
Df Sum Sq Mean Sq F value Pr(>F)
Sens 1 7906.2 7906.2 23.0210 4.119e-05 ***
Treat 2 7688.6 3844.3 11.1938 0.0002336 ***
Sens:Treat 2 1908.5 954.3 2.7786 0.0781461 .
Residuals 30 10303.0 343.4
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1$emmeans
Sens Treat emmean SE df lower.CL upper.CL
No Air 24.3 7.57 30 8.87 39.8
Yes Air 68.2 7.57 30 52.74 83.6
No Act 45.1 7.57 30 29.62 60.5
Yes Act 80.5 7.57 30 65.02 95.9
No Benz 22.2 7.57 30 6.74 37.6
Yes Benz 31.8 7.57 30 16.37 47.3
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
No Air - Yes Air -43.88 10.7 30 -4.101 0.0036
No Air - No Act -20.75 10.7 30 -1.939 0.3993
No Air - Yes Act -56.16 10.7 30 -5.249 0.0002
No Air - No Benz 2.12 10.7 30 0.199 1.0000
No Air - Yes Benz -7.51 10.7 30 -0.702 0.9803
Yes Air - No Act 23.12 10.7 30 2.161 0.2850
Yes Air - Yes Act -12.28 10.7 30 -1.148 0.8572
Yes Air - No Benz 46.00 10.7 30 4.299 0.0021
Yes Air - Yes Benz 36.37 10.7 30 3.399 0.0215
No Act - Yes Act -35.41 10.7 30 -3.309 0.0268
No Act - No Benz 22.88 10.7 30 2.138 0.2959
No Act - Yes Benz 13.24 10.7 30 1.238 0.8150
Yes Act - No Benz 58.28 10.7 30 5.447 <0.0001
Yes Act - Yes Benz 48.65 10.7 30 4.547 0.0011
[ reached 'max' / getOption("max.print") -- omitted 1 row ]
P value adjustment: tukey method for comparing a family of 6 estimates $emmeans
Sens Treat emmean SE df lower.CL upper.CL
No Air 24.3 7.57 30 8.87 39.8
Yes Air 68.2 7.57 30 52.74 83.6
No Act 45.1 7.57 30 29.62 60.5
Yes Act 80.5 7.57 30 65.02 95.9
No Benz 22.2 7.57 30 6.74 37.6
Yes Benz 31.8 7.57 30 16.37 47.3
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Yes Air - No Air 43.88 10.7 30 4.101 0.0013
No Act - No Air 20.75 10.7 30 1.939 0.2205
Yes Act - No Air 56.16 10.7 30 5.249 <0.0001
No Benz - No Air -2.12 10.7 30 -0.199 0.9958
Yes Benz - No Air 7.51 10.7 30 0.702 0.8862
P value adjustment: dunnettx method for 5 tests Otro Ejemplo El interés del estudio de Hartman-Maeir et al. fue estimar los perfiles del déficit de concienciación entre pacientes con infarto cerebral que están en rehabilitación. Estudió 35 pacientes con lesiones por infarto en el hemisferio derecho y 19 pacientes con lesiones en el hemisferio izquierdo. Además agrupó las lesiones por su tamaño como: 2 = 1-3 cm, 3 = 3-5cm, y 4 = 5 cm o más grandes. Una de las medidas importantes fue la calificación de la concienciación por su propia limitación. Las calificaciones tuvieron un rango de 8 a 24, con una calificación más alta significando mayor concienciación (REV_C08_22.csv). Pruebe la diferencia por tamaño de lesión y lado de hemisferio, sea \(\alpha = .05\).
Rev.8.22 <- read_csv("~/Dropbox/GitHub/ProbEstad/DataSets/ch08_all/REV_C08_22.csv",
show_col_types = FALSE)
boxplot(SCORES ~ factor(SIDE) + factor(SIZE), data = Rev.8.22)
Analysis of Variance Table
Response: SCORES
Df Sum Sq Mean Sq F value Pr(>F)
factor(SIZE) 2 28.690 14.3451 3.5094 0.03748 *
factor(SIDE) 1 2.578 2.5782 0.6307 0.43084
Residuals 50 204.380 4.0876
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1$emmeans
SIDE emmean SE df lower.CL upper.CL
L 11.1 0.470 50 10.11 12.0
R 10.6 0.345 50 9.91 11.3
Results are averaged over the levels of: SIZE
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
L - R 0.459 0.578 50 0.794 0.4308
Results are averaged over the levels of: SIZE Anova Table (Type II tests)
Response: SCORES
Sum Sq Df F value Pr(>F)
factor(SIZE) 29.981 2 3.7121 0.0317 *
factor(SIDE) 2.578 1 0.6384 0.4282
factor(SIZE):factor(SIDE) 10.540 2 1.3050 0.2806
Residuals 193.840 48
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1$emmeans
SIZE emmean SE df lower.CL upper.CL
2 10.0 0.429 48 9.16 10.9
3 10.5 0.560 48 9.36 11.6
4 12.1 0.535 48 11.01 13.2
Results are averaged over the levels of: SIDE
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
SIZE3 - SIZE2 0.469 0.706 48 0.664 0.7243
SIZE4 - SIZE2 2.063 0.686 48 3.009 0.0081
Results are averaged over the levels of: SIDE
P value adjustment: dunnettx method for 2 tests # emmip(Rev.8.22_lm_int, SIDE ~ SIZE)
# With rstatix
Rev.8.22_aov <- anova_test(data = Rev.8.22, SCORES ~ factor(SIZE) * factor(SIDE))
get_anova_table(Rev.8.22_aov)ANOVA Table (type II tests)
Effect DFn DFd F p p<.05 ges
1 factor(SIZE) 2 48 3.712 0.032 * 0.134
2 factor(SIDE) 1 48 0.638 0.428 0.013
3 factor(SIZE):factor(SIDE) 2 48 1.305 0.281 0.052# A tibble: 3 × 9
term .y. group1 group2 df statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 SIZE SCORES 2 3 51 -0.837 0.407 0.407 ns
2 SIZE SCORES 2 4 51 -2.65 0.0108 0.0323 *
3 SIZE SCORES 3 4 51 -1.56 0.125 0.187 ns Ejemplo 18 Los efectos de la contaminación térmica sobre Corbicula fluminea (almejas asiáticas) en tres ubicaciones geográficas diferentes fueron analizados por John Brooker (REV_C08_18). Los datos muestrales sobre la longitud, el ancho y la altura de las conchas se presentan en la siguiente tabla. Determine si existe una diferencia significativa en la longitud, altura o ancho medios (medidos en mm) de la concha en las tres ubicaciones, realizando tres análisis. ¿Qué inferencias pueden extraerse de sus resultados? ¿Cuáles son los supuestos subyacentes a dichas inferencias? ¿Cuáles son las poblaciones objetivo?
Rev.8.18 <- read_csv("~/Dropbox/GitHub/ProbEstad/DataSets/ch08_all/REV_C08_18_mod.csv",
show_col_types = FALSE)
Rev.8.18 <- Rev.8.18 %>%
mutate(Geometry = factor(Geometry))
Rev.8.18 <- Rev.8.18 %>%
mutate(site = factor(site))
ggplot(data = Rev.8.18, aes(x = site, y = Measure, fill = Geometry)) + geom_boxplot(position = position_dodge(width = 0.6),
alpha = 0.5) + geom_jitter(position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.6), aes(color = Geometry), alpha = 0.8, size = 1.4) + labs(x = "Site",
y = "Measure", fill = "Geometry") + theme_bw()
Anova Table (Type II tests)
Response: Measure
Sum Sq Df F value Pr(>F)
factor(site) 1.869 2 4.9802 0.007694 **
factor(Geometry) 286.132 2 762.5696 < 2.2e-16 ***
factor(site):factor(Geometry) 0.027 4 0.0355 0.997576
Residuals 39.961 213
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1$emmeans
site emmean SE df lower.CL upper.CL
1 5.74 0.0500 213 5.64 5.84
2 5.57 0.0521 213 5.46 5.67
3 5.54 0.0490 213 5.44 5.63
Results are averaged over the levels of: Geometry
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
site2 - site1 -0.177 0.0723 213 -2.446 0.0293
site3 - site1 -0.206 0.0700 213 -2.942 0.0071
Results are averaged over the levels of: Geometry
P value adjustment: dunnettx method for 2 tests ANOVA Modelos de medidas repetidas
Los modelos de medidas repetidas son aquellos en los que la misma
variable es medida en los mismos sujetos en dos o más ocasiones. En
donde cada ocasión puede significar una condición distinta; como
diferentes tratamientos o diferentes puntos de tiempo.
Las ventajas de usar este tipo de modelos se han mencionado y son que
permiten controlar variables relacionadas con el sujeto y con factores
relacionados con el estudio que se está desarrollado .
Ejemplo Licciardone et al. examinaron sujetos con dolor lumbar crónico inespecífico. En este estudio, 18 de los sujetos completaron un cuestionario de encuesta que evaluaba el funcionamiento físico al inicio, y después de 1, 3 y 6 meses. La Tabla (EXA_C08_S04_01) muestra los datos de estos sujetos que recibieron un tratamiento simulado que parecía ser una manipulación osteopática genuina. Valores más altos indican un mejor funcionamiento físico. El objetivo del experimento fue determinar si los sujetos reportarían mejoría con el tiempo, aun cuando el tratamiento recibido proporcionaría una mejoría mínima. Deseamos saber si existe una diferencia en los valores medios de la encuesta entre los cuatro puntos en el tiempo.
Exa8.4.1 <- read_csv("~/Dropbox/GitHub/ProbEstad/DataSets/ch08_all/EXA_C08_S04_01mod.csv",
show_col_types = FALSE)
Exa8.4.1 <- Exa8.4.1 %>%
mutate(Subject = factor(Subject), Time = factor(Time))
head(Exa8.4.1, 26)# A tibble: 26 × 3
Assessment Time Subject
<dbl> <fct> <fct>
1 80 Baseline 1
2 95 Baseline 2
3 65 Baseline 3
4 50 Baseline 4
5 60 Baseline 5
6 70 Baseline 6
7 80 Baseline 7
8 70 Baseline 8
9 80 Baseline 9
10 65 Baseline 10
# ℹ 16 more rowsggplot(data = Exa8.4.1, aes(x = Time, y = Assessment)) + geom_boxplot() + geom_jitter(width = 0.1,
alpha = 0.6, aes(color = Time)) + theme_bw()
# A tibble: 2 × 5
Time Assessment Subject is.outlier is.extreme
<fct> <dbl> <fct> <lgl> <lgl>
1 Month3 20 16 TRUE FALSE
2 Month6 25 16 TRUE FALSE # A tibble: 4 × 4
Time variable statistic p
<fct> <chr> <dbl> <dbl>
1 Baseline Assessment 0.900 0.0568
2 Month1 Assessment 0.960 0.611
3 Month3 Assessment 0.884 0.0303
4 Month6 Assessment 0.935 0.235 # # # # # # # # # # # # # # # # # # # With the R basic anova function, how to
# make a repeated measures
Exa8.4.1_mod <- aov(Assessment ~ Time + Error(Subject), data = Exa8.4.1)
summary(Exa8.4.1_mod)
Error: Subject
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 17 20238 1190
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Time 3 2396 798.6 5.501 0.00237 **
Residuals 51 7404 145.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1$emmeans
Time emmean SE df lower.CL upper.CL
Baseline 61.9 4.75 30.4 52.2 71.6
Month1 59.4 4.75 30.4 49.7 69.1
Month3 73.6 4.75 30.4 63.9 83.3
Month6 70.0 4.75 30.4 60.3 79.7
Warning: EMMs are biased unless design is perfectly balanced
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Baseline - Month1 2.50 4.02 51 0.622 0.9244
Baseline - Month3 -11.67 4.02 51 -2.905 0.0269
Baseline - Month6 -8.06 4.02 51 -2.006 0.1994
Month1 - Month3 -14.17 4.02 51 -3.527 0.0048
Month1 - Month6 -10.56 4.02 51 -2.628 0.0534
Month3 - Month6 3.61 4.02 51 0.899 0.8053
P value adjustment: tukey method for comparing a family of 4 estimates Ejemplo 2. Polyzogopoulou et al. reportan los efectos de la cirugía bariátrica sobre los niveles de glucosa en ayuno (mmol/L) en 12 sujetos obesos con diabetestipo 2 en cuatro momentos: antes de la operación, a los 3 meses, 6 meses y 12 meses. ¿Podemos concluir, tras eliminar los efectos atribuibles a los sujetos, que los niveles de glucosa en ayuno difieren a lo largo del tiempo después de la cirugía? Considere \(\alpha = 0.05\).
Exr8.32 <- read_csv("~/Dropbox/GitHub/ProbEstad/DataSets/ch08_all/REV_C08_32.csv",
show_col_types = FALSE)
Exr8.32# A tibble: 48 × 3
SUBJ GLUCOSE TIME
<dbl> <dbl> <dbl>
1 1 108 0
2 2 96.7 0
3 3 77 0
4 4 92 0
5 5 97 0
6 6 94 0
7 7 76 0
8 8 100 0
9 9 82 0
10 10 104. 0
# ℹ 38 more rowsExr8.32 <- Exr8.32 %>%
mutate(SUBJ = factor(SUBJ), TIME = factor(TIME))
ggplot(data = Exr8.32, aes(x = TIME, y = GLUCOSE)) + geom_boxplot() + geom_jitter(width = 0.1,
alpha = 0.6, aes(color = TIME)) + theme_bw()
# A tibble: 1 × 5
TIME SUBJ GLUCOSE is.outlier is.extreme
<fct> <fct> <dbl> <lgl> <lgl>
1 3 9 282 TRUE FALSE # A tibble: 4 × 4
TIME variable statistic p
<fct> <chr> <dbl> <dbl>
1 0 GLUCOSE 0.940 0.497
2 3 GLUCOSE 0.902 0.169
3 6 GLUCOSE 0.951 0.652
4 12 GLUCOSE 0.925 0.334
Error: SUBJ
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 11 10517 956.1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
TIME 3 65135 21712 40.22 3.98e-11 ***
Residuals 33 17816 540
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1$emmeans
TIME emmean SE df lower.CL upper.CL
0 90.5 7.33 40.8 75.7 105.3
3 170.8 7.33 40.8 156.0 185.6
6 84.7 7.33 40.8 69.9 99.5
12 82.8 7.33 40.8 68.0 97.6
Warning: EMMs are biased unless design is perfectly balanced
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
TIME0 - TIME3 -80.32 9.49 33 -8.467 <0.0001
TIME0 - TIME6 5.83 9.49 33 0.614 0.9269
TIME0 - TIME12 7.68 9.49 33 0.810 0.8493
TIME3 - TIME6 86.14 9.49 33 9.081 <0.0001
TIME3 - TIME12 88.00 9.49 33 9.277 <0.0001
TIME6 - TIME12 1.86 9.49 33 0.196 0.9973
P value adjustment: tukey method for comparing a family of 4 estimates $emmeans
TIME emmean SE df lower.CL upper.CL
0 90.5 7.33 40.8 75.7 105.3
3 170.8 7.33 40.8 156.0 185.6
6 84.7 7.33 40.8 69.9 99.5
12 82.8 7.33 40.8 68.0 97.6
Warning: EMMs are biased unless design is perfectly balanced
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
TIME3 - TIME0 80.32 9.49 33 8.467 <0.0001
TIME6 - TIME0 -5.83 9.49 33 -0.614 0.8439
TIME12 - TIME0 -7.68 9.49 33 -0.810 0.7359
P value adjustment: dunnettx method for 3 tests Hay un par de librerias de R - rstatix y ggpubr - que redefinen varias de las puebas de “base R” y de ggplot para estructuras tibble. En esta sección usamos algunas de estas funciones que en ocaciones hacen más fácil la definición del problema y el uso de los datos en tibble.
# ggplot used with ggpubr package
bxp <- ggboxplot(data = Exa8.4.1, x = "Time", y = "Assessment", add = "jitter", color = "Time")
bxp
# ANOVA test from the rstatix package
Exa8.4.1_aov <- anova_test(data = Exa8.4.1, dv = Assessment, wid = Subject, within = Time)
get_anova_table(Exa8.4.1_aov)ANOVA Table (type III tests)
Effect DFn DFd F p p<.05 ges
1 Time 2.22 37.68 5.501 0.006 * 0.08# A tibble: 6 × 9
term .y. group1 group2 df statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 Time Assessment Baseline Month1 68 0.372 0.711 1 ns
2 Time Assessment Baseline Month3 68 -1.74 0.0871 0.523 ns
3 Time Assessment Baseline Month6 68 -1.20 0.235 1 ns
4 Time Assessment Month1 Month3 68 -2.11 0.0387 0.232 ns
5 Time Assessment Month1 Month6 68 -1.57 0.121 0.725 ns
6 Time Assessment Month3 Month6 68 0.537 0.593 1 ns # Pairwise comparisons
pwc <- Exa8.4.1 %>%
pairwise_t_test(Assessment ~ Time, paired = TRUE, p.adjust.method = "bonferroni")
pwc# A tibble: 6 × 10
.y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
* <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
1 Assessment Baseli… Month1 18 18 0.658 17 0.519 1 ns
2 Assessment Baseli… Month3 18 18 -2.72 17 0.015 0.088 ns
3 Assessment Baseli… Month6 18 18 -1.75 17 0.099 0.592 ns
4 Assessment Month1 Month3 18 18 -3.86 17 0.001 0.007 **
5 Assessment Month1 Month6 18 18 -2.22 17 0.04 0.24 ns
6 Assessment Month3 Month6 18 18 1.40 17 0.18 1 ns