Bibliotecas

library(tidyverse)
library(effectsize)
library(readxl)
library(broom)
library(emmeans)
library(nlme)

Datos

datos <- read_excel("experimento_ratones.xlsx") 
datos

Tamaño del efecto

Modelo lineal

Modelo

modelo6 <- aov(GST ~ Treatment * Strain + Block, data = datos)
modelo6 %>% tidy()

Medias estimadas

  • Tratamiento:-
modelo6 %>% 
  emmeans(specs = "Treatment") %>% 
  as.data.frame()
  • Raza:
modelo6 %>% 
  emmeans(specs = "Strain") %>% 
  as.data.frame()
  • Interacción raza-tratamiento:
modelo6 %>% 
  emmeans(specs = pairwise ~ Treatment | Strain) %>% 
  as.data.frame()

Tamaño del efecto

modelo6 %>% 
  eta_squared()

Residuales

par(mfrow = c(2, 2))
plot(modelo6)

Modelo lineal mixto

Modelo

modelo_mixto1 <- lme(GST ~ Treatment * Strain,
                     random = ~ 1 | Block,
                     data = datos)

modelo_mixto1 %>% anova()

Diferencia de medias

  • Interacción raza-tratamiento:
modelo_mixto1 %>% 
  emmeans(specs = pairwise ~ Treatment | Strain) %>% 
  pairs() %>% 
  tidy()

Tamaño del efecto

modelo_mixto1 %>% 
  eta_squared()

Residuales

  • Gráfico de normalidad:
residuales <- residuals(modelo_mixto1)
reales <- datos$GST
ajustados <- fitted(modelo_mixto1)

ggpubr::ggqqplot(residuales)

shapiro.test(residuales)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.93226, p-value = 0.2647
  • Homocedasticidad:
plot(x = ajustados, y = residuales)
abline(h = 0, lty = 2, col = "red")

Comparación de modelos

AIC(modelo6, modelo_mixto1)
BIC(modelo6, modelo_mixto1)

Medidas repetidas

Datos de ejemplo

df_peso <- ChickWeight 
df_peso %>% head()

Gráfico exploratorio

df_peso %>%
  ggplot(aes(x = Time, y = weight, color = Diet)) +
  geom_point() +
  theme_bw() +
  facet_wrap( ~ Chick)

Modelo lineal 1

modelo_pesos1 <- aov(weight ~ Time + Diet, data = df_peso)
modelo_pesos1 %>% summary()
##              Df  Sum Sq Mean Sq F value Pr(>F)    
## Time          1 2042344 2042344 1576.46 <2e-16 ***
## Diet          3  129876   43292   33.42 <2e-16 ***
## Residuals   573  742336    1296                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelo lineal 2

modelo_pesos2 <- aov(weight ~ Time * Diet, data = df_peso)
modelo_pesos2 %>% summary()
##              Df  Sum Sq Mean Sq F value   Pr(>F)    
## Time          1 2042344 2042344 1759.76  < 2e-16 ***
## Diet          3  129876   43292   37.30  < 2e-16 ***
## Time:Diet     3   80804   26935   23.21 3.47e-14 ***
## Residuals   570  661532    1161                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelos mixtos

modelo_pesos3 <- lme(weight ~ Time + Diet, data = df_peso, random = ~ 1 | Chick)
modelo_pesos4 <- lme(weight ~ Time + Diet, data = df_peso, random = ~ Time | Chick)
modelo_pesos5 <- lme(weight ~ Time * Diet, data = df_peso, random = ~ Time | Chick)

modelo_pesos3 %>% anova()
modelo_pesos4 %>% anova()
modelo_pesos5 %>% anova()

Comparación de modelos

AIC(modelo_pesos1, modelo_pesos2, modelo_pesos3, modelo_pesos4, modelo_pesos5)
BIC(modelo_pesos1, modelo_pesos2, modelo_pesos3, modelo_pesos4, modelo_pesos5)

Medias estimadas

  • Interacción Time-Dieta - Gráfico:
emmip(modelo_pesos5, Diet ~ Time , cov.reduce = range)

Tamaño del efecto

modelo_pesos5 %>% 
  eta_squared()

Residuales

  • Gráfico de normalidad:
residuales <- residuals(modelo_pesos5)
reales <- df_peso$weight
ajustados <- fitted(modelo_pesos5)

ggpubr::ggqqplot(residuales)

shapiro.test(residuales)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.9873, p-value = 6.39e-05
  • Homocedasticidad:
plot(x = ajustados, y = residuales)
abline(h = 0, lty = 2, col = "red")