Experimento con 1 factor

Datos

library(tidyverse)
library(readxl)
library(janitor)
datos_cafe <- read_excel("datos_café.xlsx") %>% 
  clean_names() %>% 
  mutate(mycorrhizal_colonization = as.numeric(mycorrhizal_colonization))
datos_cafe

## # A tibble: 81 × 4
##    system period localitie mycorrhizal_colonization
##    <chr>  <chr>  <chr>                        <dbl>
##  1 Forest 1st    Z1                            60  
##  2 Forest 1st    Z1                            50  
##  3 Forest 1st    Z1                            31.3
##  4 Forest 1st    Z2                            10.3
##  5 Forest 1st    Z2                            20  
##  6 Forest 1st    Z2                            28.6
##  7 Forest 1st    Z3                            49  
##  8 Forest 1st    Z3                            42  
##  9 Forest 1st    Z3                            50  
## 10 Forest 2nd    Z1                             9  
## # … with 71 more rows

Resumen descriptivo básico

datos_cafe %>% 
  group_by(system) %>% 
  summarise(promedio = mean(mycorrhizal_colonization),
            desviacion = sd(mycorrhizal_colonization),
            N = n())

## # A tibble: 3 × 4
##   system         promedio desviacion     N
##   <chr>             <dbl>      <dbl> <int>
## 1 Agroecological     25.0      15.0     27
## 2 Conventional       36.3       9.97    27
## 3 Forest             25.5      16.9     27

Boxplot

datos_cafe %>% 
  ggplot(mapping = aes(x = system, y = mycorrhizal_colonization)) +
  geom_boxplot()

Modelo estadístico

\[y_{ij} = \mu + \beta_j + \epsilon_{ij}\]

Donde:

\(y_{ij}\): variable respuesta. En este caso es la colonización micorrizal.
\(\mu\): media general para la variable respuesta.
\(\beta_j\): efecto del \(j-ésimo\) sistema sobre la variable colonización micorrizal. Donde \(j = 1, 2, 3\)
\(\epsilon_{ij}\): error experimental

Hipótesis

\[H_0: \mu_{Conventional} = \mu_{Agroecological} = \mu_{Forest} \\ H_1: \mu_i \neq \mu_j\]

Nivel de significancia

En este caso vamos a usar el 5%.

Anova

En R podemos usar la función “aov()” para ejecutar análisis de varianza.
También podemos utilizar la función “lm()” para ejecutar análisis de varianza.

modelo_cafe <- aov(formula = mycorrhizal_colonization ~ system,
                   data = datos_cafe)

# Resumen del modelo
summary(modelo_cafe)

##             Df Sum Sq Mean Sq F value Pr(>F)   
## system       2   2211  1105.7    5.43 0.0062 **
## Residuals   78  15884   203.6                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusión: como el valor p (0.0062) es menor que el nivel de significancia (0.05), existe evidencia para rechazar la hipótesis nula, es decir, que al menos un par de sistemas difieren estadísticamente para la colonización micorrizal.

Análisis de residuales

Podemos obtener los residuales del modelo con la función “residuals()”:

residuales_cafe <- residuals(modelo_cafe)
residuales_cafe

##           1           2           3           4           5           6 
##  34.5333333  24.5333333   5.8333333 -15.1666667  -5.4666667   3.1333333 
##           7           8           9          10          11          12 
##  23.5333333  16.5333333  24.5333333 -16.4666667 -14.4666667 -11.4666667 
##          13          14          15          16          17          18 
## -13.1666667   1.8333333 -18.1666667 -11.1666667 -16.4666667  -9.7666667 
##          19          20          21          22          23          24 
##   4.2333333 -18.7666667 -16.4666667 -13.7666667  -7.7666667  26.5333333 
##          25          26          27          28          29          30 
##   3.2333333  -6.1666667  26.2333333  39.2629630  42.2629630  -8.4370370 
##          31          32          33          34          35          36 
##  -3.7370370  -2.7370370   6.6629630   4.2629630   0.2629630 -13.7370370 
##          37          38          39          40          41          42 
## -16.0370370 -10.7370370 -19.7370370 -16.7370370 -15.0370370  -0.3370370 
##          43          44          45          46          47          48 
## -12.3370370   3.2629630  -1.3370370   3.9629630   6.9629630 -10.3370370 
##          49          50          51          52          53          54 
##  -0.3370370 -10.7370370  -0.7370370  10.9629630  13.2629630  11.9629630 
##          55          56          57          58          59          60 
## -11.0296296   9.6703704  26.9703704  -6.3296296   8.6703704  -5.0296296 
##          61          62          63          64          65          66 
## -13.3296296   7.6703704  -7.6296296  -7.3296296   6.9703704  -9.3296296 
##          67          68          69          70          71          72 
## -13.6296296  -5.6296296   0.9703704  16.6703704  -7.0296296  -5.6296296 
##          73          74          75          76          77          78 
##  -9.3296296   2.9703704   8.6703704  11.9703704   0.6703704  -9.3296296 
##          79          80          81 
##   5.3703704   0.3703704   2.9703704

Normalidad

library(ggpubr)
ggqqplot(residuales_cafe)

Prueba de shapiro wilk: la normalidad no se cumple, sin embargo, los desvíos que se observan en el gŕafico parecen ser pequeños.

shapiro.test(residuales_cafe)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_cafe
## W = 0.92287, p-value = 0.0001182

Homocedasticidad

La homocedasticidad sí se cumple.

library(car)
leveneTest(residuales_cafe ~ datos_cafe$system)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  2  2.3892 0.09838 .
##       78                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Gráfico de ajustados vs residuales:

ajustados_cafe <- fitted(modelo_cafe)

plot(x = ajustados_cafe, y = residuales_cafe)

Gráficos de residuales por defecto con la función “plot()”:

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

Estimaciones

Promedios

Con la función “model.tables()” podemos obtener los promedios estimados:

model.tables(x = modelo_cafe,
             type = "means",
             se = TRUE)

## Tables of means
## Grand mean
##          
## 28.94444 
## 
##  system 
## system
## Agroecological   Conventional         Forest 
##          25.04          36.33          25.47 
## 
## Standard errors for differences of means
##         system
##          3.884
## replic.     27

Efectos

Los efectos (positivos o negativos) pueden ser obtenidos con la misma función “model.tables()”

model.tables(x = modelo_cafe,
             type = "effects",
             se = TRUE)

## Tables of effects
## 
##  system 
## system
## Agroecological   Conventional         Forest 
##         -3.907          7.385         -3.478 
## 
## Standard errors of effects
##         system
##          2.746
## replic.     27

Promedios e IC

En este caso vamos a usar la biblioteca ggeffects:

library(ggeffects)
ggeffect(model = modelo_cafe)

## $system
## # Predicted values of mycorrhizal_colonization
## 
## system         | Predicted |         95% CI
## -------------------------------------------
## Agroecological |     25.04 | [19.57, 30.50]
## Conventional   |     36.33 | [30.86, 41.80]
## Forest         |     25.47 | [20.00, 30.93]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "modelo_cafe"

Comparaciones múltiples

Podemos utilizar métodos con y sin ajuste del valor p.

Sin ajuste

pairwise.t.test(x = datos_cafe$mycorrhizal_colonization,
                g = datos_cafe$system,
                p.adjust.method = "none", 
                alternative = "two.sided")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  datos_cafe$mycorrhizal_colonization and datos_cafe$system 
## 
##              Agroecological Conventional
## Conventional 0.0047         -           
## Forest       0.9122         0.0065      
## 
## P value adjustment method: none

Con ajuste

Bonferroni:

pairwise.t.test(x = datos_cafe$mycorrhizal_colonization,
                g = datos_cafe$system,
                p.adjust.method = "bonferroni", 
                alternative = "two.sided")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  datos_cafe$mycorrhizal_colonization and datos_cafe$system 
## 
##              Agroecological Conventional
## Conventional 0.014          -           
## Forest       1.000          0.019       
## 
## P value adjustment method: bonferroni

Holm:

pairwise.t.test(x = datos_cafe$mycorrhizal_colonization,
                g = datos_cafe$system,
                p.adjust.method = "holm", 
                alternative = "two.sided")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  datos_cafe$mycorrhizal_colonization and datos_cafe$system 
## 
##              Agroecological Conventional
## Conventional 0.014          -           
## Forest       0.912          0.014       
## 
## P value adjustment method: holm

Tukey

TukeyHSD(x = modelo_cafe)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mycorrhizal_colonization ~ system, data = datos_cafe)
## 
## $system
##                                    diff        lwr       upr     p adj
## Conventional-Agroecological  11.2925926   2.012968 20.572217 0.0130187
## Forest-Agroecological         0.4296296  -8.849995  9.709254 0.9932772
## Forest-Conventional         -10.8629630 -20.142587 -1.583339 0.0176515

Reporte de resultados

Gráficos

Medias estimadas e intervalos de confianza:

g1 <- ggeffect(model = modelo_cafe)
plot(g1)

## $system

Gráfico de comparación de medias con Tukey:

library(broom)
g2 <- TukeyHSD(x = modelo_cafe) %>% tidy()
g2 %>% 
  ggplot(mapping = aes(x = contrast, y = estimate)) +
  geom_point() +
  geom_errorbar(mapping = aes(ymin = conf.low, ymax = conf.high),
                width = 0.2) +
  geom_hline(yintercept = 0, lty = 2, color = "red") +
  coord_flip()

Tamaño del efecto

En este caso el tamaño del efecto es moderado o medio.

library(effectsize)
eta_squared(model = modelo_cafe)

## # Effect Size for ANOVA
## 
## Parameter | Eta2 |       90% CI
## -------------------------------
## system    | 0.12 | [0.02, 0.23]

DBCA

DBCA: Diseño en Bloques Completamente Aleatorizados

Datos

fertilidad <- read_excel("fertilidad.xlsx") %>%
  mutate(
    tratamiento = factor(
      tratamiento,
      levels = c("Control", "CaNO32", "CONH22",
                 "NaNO3", "NH42SO4", "NH4NO3")
    ),
    nombre = factor(
      nombre,
      levels = c(
        "Sin nitrógeno",
        "Nitrato de calcio",
        "Urea",
        "Nitrato de sodio",
        "Sulfato de amonio",
        "Nitrato de amonio"
      )
    ),
    tipo_suelo = factor(tipo_suelo,
                        levels = c("I", "II", "III", "IV"))
  )
fertilidad

## # A tibble: 24 × 4
##    tratamiento tipo_suelo produccion nombre           
##    <fct>       <fct>           <dbl> <fct>            
##  1 NH42SO4     I                32.1 Sulfato de amonio
##  2 NH42SO4     II               35.6 Sulfato de amonio
##  3 NH42SO4     III              41.9 Sulfato de amonio
##  4 NH42SO4     IV               35.4 Sulfato de amonio
##  5 NH4NO3      I                30.1 Nitrato de amonio
##  6 NH4NO3      II               31.5 Nitrato de amonio
##  7 NH4NO3      III              37.1 Nitrato de amonio
##  8 NH4NO3      IV               30.8 Nitrato de amonio
##  9 CONH22      I                25.4 Urea             
## 10 CONH22      II               27.1 Urea             
## # … with 14 more rows

Estadísticos descriptivos

Por Tratamiento

fertilidad %>% 
  group_by(nombre) %>% 
  summarise(promedio = mean(produccion),
            desviacion = sd(produccion),
            N = n())

## # A tibble: 6 × 4
##   nombre            promedio desviacion     N
##   <fct>                <dbl>      <dbl> <int>
## 1 Sin nitrógeno         25.4       1.69     4
## 2 Nitrato de calcio     31.0       4.93     4
## 3 Urea                  29.4       3.81     4
## 4 Nitrato de sodio      30.7       3.28     4
## 5 Sulfato de amonio     36.2       4.09     4
## 6 Nitrato de amonio     32.4       3.20     4

Por Bloque

fertilidad %>% 
  group_by(tipo_suelo) %>% 
  summarise(promedio = mean(produccion),
            desviacion = sd(produccion),
            N = n())

## # A tibble: 4 × 4
##   tipo_suelo promedio desviacion     N
##   <fct>         <dbl>      <dbl> <int>
## 1 I              26.8       3.51     6
## 2 II             30.5       3.94     6
## 3 III            34.8       4.98     6
## 4 IV             31.2       2.78     6

Gráfico

fertilidad %>% 
  ggplot(mapping = aes(x = nombre, y = produccion, color = tipo_suelo)) +
  geom_point() +
  geom_line(aes(group = tipo_suelo)) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Modelo

\[ y_{ij} = \mu + \alpha_i + \beta_j + \epsilon_{ij}, \qquad i = 1, 2, \cdots, 6, \quad j = 1, 2, \cdots, 4 \\ \epsilon_{ij} \sim \mathcal{N}(0,\sigma^2) \quad i.i.d. \]

Hipótesis

Principal

\[ H_0: \alpha_1 = \alpha_2 = \alpha_3 = \alpha_4 = \alpha_5 = \alpha_6 = 0 \\ H_A: \textrm{Alguna } \alpha_i \textrm{ diferente de } 0 \]

Verificación

\[ H_0: \beta_1 = \beta_2 = \beta_3 = \beta_4 = 0 \\ H_A: \textrm{Alguna } \beta_i \textrm{ diferente de } 0 \]

Nivel de significancia

En este caso vamos a usar el 5% (0.05)

Anova

Conclusiones:
- Como el valor p (1.10e-05) para la fuente de variación “nombre” es menor que el nivel de significancia (0.05), existe evidencia para rechazar la hipótesis nula, es decir, que al menos un par de fuentes de nitrógeno discrepan estadísticamente.
- Como el valor p (1.22e-05 ) para la fuente de variación “tipo_suelo” es menor que el nivel de significancia (0.05), existe evidencia para rechazar la hipótesis nula (hipótesis de verificación), es decir, que al menos un par de tipos de suelo discrepan estadísticamente.

modelo1 <- aov(produccion ~ nombre + tipo_suelo,
               data = fertilidad)
summary(modelo1)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## nombre       5 256.15   51.23   16.85 1.10e-05 ***
## tipo_suelo   3 192.75   64.25   21.13 1.22e-05 ***
## Residuals   15  45.62    3.04                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Análisis de residuales

residuales <- residuals(modelo1)

Normalidad

Conclusión: como el valor p (0.6013) es mayor que el nivel de significancia, no existe evidencia para rechazar que los residuales se distribuyen de forma normal.

shapiro.test(residuales)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.96731, p-value = 0.6013

ggqqplot(residuales)

Homocedasticidad

Tratamiento: sí se cumple la homocedasticidad.

bartlett.test(residuales ~ fertilidad$nombre)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuales by fertilidad$nombre
## Bartlett's K-squared = 2.715, df = 5, p-value = 0.7438

Bloque: sí se cumple la homocedasticidad.

bartlett.test(residuales ~ fertilidad$tipo_suelo)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuales by fertilidad$tipo_suelo
## Bartlett's K-squared = 0.4122, df = 3, p-value = 0.9377

Estimaciones

Promedios

model.tables(x = modelo1,
             type = "means",
             se = TRUE)

## Tables of means
## Grand mean
##          
## 30.84167 
## 
##  nombre 
## nombre
##     Sin nitrógeno Nitrato de calcio              Urea  Nitrato de sodio 
##             25.35             31.03             29.35             30.70 
## Sulfato de amonio Nitrato de amonio 
##             36.25             32.38 
## 
##  tipo_suelo 
## tipo_suelo
##     I    II   III    IV 
## 26.83 30.50 34.82 31.22 
## 
## Standard errors for differences of means
##         nombre tipo_suelo
##          1.233      1.007
## replic.      4          6

Efectos

model.tables(x = modelo1,
             type = "effects",
             se = TRUE)

## Tables of effects
## 
##  nombre 
## nombre
##     Sin nitrógeno Nitrato de calcio              Urea  Nitrato de sodio 
##            -5.492             0.183            -1.492            -0.142 
## Sulfato de amonio Nitrato de amonio 
##             5.408             1.533 
## 
##  tipo_suelo 
## tipo_suelo
##      I     II    III     IV 
## -4.008 -0.342  3.975  0.375 
## 
## Standard errors of effects
##         nombre tipo_suelo
##         0.8719     0.7119
## replic.      4          6

Promedios e IC

Conclusión:
- Dada la construcción matemática podemos afirmar con un nivel de confianza del 95%, que en esta muestra la media de producción para el tratamiento “sin nitrógeno” fue de 25.35 (estimación puntual), sin embargo, dada la variabilidad, podría esperar que la verdadera media de la producción para el tratamiento “sin nitrógeno” esté entre 23.49 y 27.21.

ggeffect(model = modelo1)

## $nombre
## # Predicted values of produccion
## 
## nombre            | Predicted |         95% CI
## ----------------------------------------------
## Sin nitrógeno     |     25.35 | [23.49, 27.21]
## Nitrato de calcio |     31.03 | [29.17, 32.88]
## Urea              |     29.35 | [27.49, 31.21]
## Nitrato de sodio  |     30.70 | [28.84, 32.56]
## Sulfato de amonio |     36.25 | [34.39, 38.11]
## Nitrato de amonio |     32.38 | [30.52, 34.23]
## 
## $tipo_suelo
## # Predicted values of produccion
## 
## tipo_suelo | Predicted |         95% CI
## ---------------------------------------
## I          |     26.83 | [25.32, 28.35]
## II         |     30.50 | [28.98, 32.02]
## III        |     34.82 | [33.30, 36.33]
## IV         |     31.22 | [29.70, 32.73]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "modelo1"

Comparaciones múltiples

Tukey

Nota: con “which” podemos establecer el nombre de la fuente de variación sobre la cual queremos aplicar las comparaciones pareadas.

TukeyHSD(modelo1, which = "nombre") %>% tidy()

## # A tibble: 15 × 7
##    term   contrast            null.value estimate conf.low conf.high adj.p.value
##    <chr>  <chr>                    <dbl>    <dbl>    <dbl>     <dbl>       <dbl>
##  1 nombre Nitrato de calcio-…          0    5.68   1.67        9.68   0.00381   
##  2 nombre Urea-Sin nitrógeno           0    4     -0.00633     8.01   0.0505    
##  3 nombre Nitrato de sodio-S…          0    5.35   1.34        9.36   0.00631   
##  4 nombre Sulfato de amonio-…          0   10.9    6.89       14.9    0.00000312
##  5 nombre Nitrato de amonio-…          0    7.03   3.02       11.0    0.000498  
##  6 nombre Urea-Nitrato de ca…          0   -1.68  -5.68        2.33   0.750     
##  7 nombre Nitrato de sodio-N…          0   -0.325 -4.33        3.68   1.00      
##  8 nombre Sulfato de amonio-…          0    5.22   1.22        9.23   0.00766   
##  9 nombre Nitrato de amonio-…          0    1.35  -2.66        5.36   0.876     
## 10 nombre Nitrato de sodio-U…          0    1.35  -2.66        5.36   0.876     
## 11 nombre Sulfato de amonio-…          0    6.90   2.89       10.9    0.000598  
## 12 nombre Nitrato de amonio-…          0    3.03  -0.981       7.03   0.200     
## 13 nombre Sulfato de amonio-…          0    5.55   1.54        9.56   0.00463   
## 14 nombre Nitrato de amonio-…          0    1.68  -2.33        5.68   0.750     
## 15 nombre Nitrato de amonio-…          0   -3.87  -7.88        0.131  0.0608

Gráfico de diferencia estimada con intervalo de confianzacdel 95%:

TukeyHSD(modelo1, which = "nombre") %>%
  tidy() %>% 
  ggplot(mapping = aes(x = contrast, y = estimate)) +
  geom_point() +
  geom_errorbar(mapping = aes(ymin = conf.low, ymax = conf.high),
                width = 0.2) +
  geom_hline(yintercept = 0, lty = 2, color = "red") +
  coord_flip()

Dunnet

library(multcomp)
PruebaDunnet <- glht(modelo1, linfct = mcp(nombre = "Dunnett"))
summary(PruebaDunnet)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Dunnett Contrasts
## 
## 
## Fit: aov(formula = produccion ~ nombre + tipo_suelo, data = fertilidad)
## 
## Linear Hypotheses:
##                                        Estimate Std. Error t value Pr(>|t|)    
## Nitrato de calcio - Sin nitrógeno == 0    5.675      1.233   4.602  0.00150 ** 
## Urea - Sin nitrógeno == 0                 4.000      1.233   3.244  0.02168 *  
## Nitrato de sodio - Sin nitrógeno == 0     5.350      1.233   4.339  0.00261 ** 
## Sulfato de amonio - Sin nitrógeno == 0   10.900      1.233   8.839  < 0.001 ***
## Nitrato de amonio - Sin nitrógeno == 0    7.025      1.233   5.697  < 0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Tamaño del efecto

Nota: cuando tenemos más de un factor, el eta cuadrado y el índice F de Cohen, son interpretados de manera parcial.

eta_squared(modelo1)

## # Effect Size for ANOVA (Type I)
## 
## Parameter  | Eta2 (partial) |       90% CI
## ------------------------------------------
## nombre     |           0.85 | [0.67, 0.91]
## tipo_suelo |           0.81 | [0.61, 0.88]

Consecuencias de no bloquear

Modelo

modelo2 <- aov(produccion ~ nombre, data = fertilidad)
summary(modelo2)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## nombre       5  256.1   51.23   3.869 0.0148 *
## Residuals   18  238.4   13.24                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tukey

TukeyHSD(modelo2) %>% tidy()

## # A tibble: 15 × 7
##    term   contrast            null.value estimate conf.low conf.high adj.p.value
##    <chr>  <chr>                    <dbl>    <dbl>    <dbl>     <dbl>       <dbl>
##  1 nombre Nitrato de calcio-…          0    5.68     -2.50     13.9      0.283  
##  2 nombre Urea-Sin nitrógeno           0    4        -4.18     12.2      0.636  
##  3 nombre Nitrato de sodio-S…          0    5.35     -2.83     13.5      0.340  
##  4 nombre Sulfato de amonio-…          0   10.9       2.72     19.1      0.00560
##  5 nombre Nitrato de amonio-…          0    7.03     -1.15     15.2      0.117  
##  6 nombre Urea-Nitrato de ca…          0   -1.68     -9.85      6.50     0.985  
##  7 nombre Nitrato de sodio-N…          0   -0.325    -8.50      7.85     1.00   
##  8 nombre Sulfato de amonio-…          0    5.22     -2.95     13.4      0.364  
##  9 nombre Nitrato de amonio-…          0    1.35     -6.83      9.53     0.994  
## 10 nombre Nitrato de sodio-U…          0    1.35     -6.83      9.53     0.994  
## 11 nombre Sulfato de amonio-…          0    6.90     -1.28     15.1      0.128  
## 12 nombre Nitrato de amonio-…          0    3.03     -5.15     11.2      0.843  
## 13 nombre Sulfato de amonio-…          0    5.55     -2.63     13.7      0.304  
## 14 nombre Nitrato de amonio-…          0    1.68     -6.50      9.85     0.985  
## 15 nombre Nitrato de amonio-…          0   -3.87    -12.1       4.30     0.665

Diseño factorial 2 vías

Datos

naranja <- read_csv2("experimento_vitamC.csv") %>% 
  mutate(dias = factor(dias)) 
naranja

## # A tibble: 36 × 3
##    dias  tipoJugo vitamC
##    <fct> <chr>     <dbl>
##  1 0     B          52.6
##  2 0     B          54.2
##  3 0     B          49.8
##  4 0     B          46.5
##  5 0     C          56  
##  6 0     C          48  
##  7 0     C          49.6
##  8 0     C          48.4
##  9 0     A          52.5
## 10 0     A          52  
## # … with 26 more rows

Estadísticos descriptivos

Tipo de jugo

naranja %>% 
  group_by(tipoJugo) %>% 
  summarise(promedio = mean(vitamC),
            desviacion = sd(vitamC),
            N = n())

## # A tibble: 3 × 4
##   tipoJugo promedio desviacion     N
##   <chr>       <dbl>      <dbl> <int>
## 1 A            49.0       3.08    12
## 2 B            48.1       4.36    12
## 3 C            46.6       4.10    12

Días

naranja %>% 
  group_by(dias) %>% 
  summarise(promedio = mean(vitamC),
            desviacion = sd(vitamC),
            N = n())

## # A tibble: 3 × 4
##   dias  promedio desviacion     N
##   <fct>    <dbl>      <dbl> <int>
## 1 0         51.2       2.81    12
## 2 3         47.2       3.27    12
## 3 7         45.2       3.01    12

Días y Tipo de jugo

naranja %>% 
  group_by(dias, tipoJugo) %>% 
  summarise(promedio = mean(vitamC),
            desviacion = sd(vitamC),
            N = n())

## # A tibble: 9 × 5
## # Groups:   dias [3]
##   dias  tipoJugo promedio desviacion     N
##   <fct> <chr>       <dbl>      <dbl> <int>
## 1 0     A            52.5      0.806     4
## 2 0     B            50.8      3.38      4
## 3 0     C            50.5      3.73      4
## 4 3     A            48.2      1.07      4
## 5 3     B            48.6      4.31      4
## 6 3     C            44.8      2.77      4
## 7 7     A            46.2      2.32      4
## 8 7     B            44.9      3.98      4
## 9 7     C            44.6      3.17      4

Gráfico 1

naranja %>% 
  ggplot(mapping = aes(x = dias, y = vitamC, color = tipoJugo, fill = tipoJugo)) +
  geom_boxplot(alpha = 0.5)

Gráfico 2

naranja %>% 
  group_by(dias, tipoJugo) %>% 
  summarise(promedio = mean(vitamC)) %>% 
  ungroup() %>% 
  ggplot(mapping = aes(x = dias, y = promedio, color = tipoJugo,
                       group = tipoJugo)) +
  geom_point(size = 2) +
  geom_line()

Modelo

\[ y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha \beta)_{ij} + \epsilon_{ijk} \\ i = 1, 2, 3. \quad j = 1, 2, 3. \quad \\ \epsilon_{ijk} \sim \mathcal{N}(0, \sigma^2) i.i.d. \]

Hipótesis

Tipo de jugo

\[ H_0: \alpha_1 = \alpha_2 = \alpha_3 = 0 \\ H_1: \textrm{Alguna } \alpha \textrm{ diferente de 0} \]

Días

\[ H_0: \beta_1 = \beta_2 = \beta_3 = 0 \\ H_1: \textrm{Algún } \beta \textrm{ diferente de 0} \]

Interacción

\[ H_0: (\alpha\beta)_{11} = \cdots = (\alpha\beta)_{33} = 0 \\ H_1: \textrm{Algún } (\alpha\beta)_{ij} \textrm{ diferente de 0} \]

Nivel de significancia

En este caso usamos un valor de 1%.

Anova

Modelo 1

m1 <- aov(vitamC ~ dias + tipoJugo + dias:tipoJugo, data = naranja)
summary(m1)

##               Df Sum Sq Mean Sq F value   Pr(>F)    
## dias           2 226.68  113.34  12.041 0.000183 ***
## tipoJugo       2  32.96   16.48   1.751 0.192749    
## dias:tipoJugo  4  17.30    4.33   0.460 0.764691    
## Residuals     27 254.14    9.41                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelo 2

m2 <- aov(vitamC ~ dias + tipoJugo, data = naranja)
summary(m2)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## dias         2 226.68  113.34  12.944 8.19e-05 ***
## tipoJugo     2  32.96   16.48   1.882    0.169    
## Residuals   31 271.44    8.76                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Análisis de residuales

residuales <- residuals(m2)

Normalidad

shapiro.test(residuales)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.97457, p-value = 0.5627

Homocedasticidad

Tipo de jugo:

bartlett.test(naranja$vitamC ~ naranja$tipoJugo)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  naranja$vitamC by naranja$tipoJugo
## Bartlett's K-squared = 1.3539, df = 2, p-value = 0.5082

Bloque:

bartlett.test(naranja$vitamC ~ naranja$dias)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  naranja$vitamC by naranja$dias
## Bartlett's K-squared = 0.24152, df = 2, p-value = 0.8862

Estimaciones

Promedios

model.tables(x = m2,
             type = "means",
             se = TRUE)

## Tables of means
## Grand mean
##          
## 47.89444 
## 
##  dias 
## dias
##     0     3     7 
## 51.25 47.22 45.22 
## 
##  tipoJugo 
## tipoJugo
##     A     B     C 
## 48.95 48.10 46.63 
## 
## Standard errors for differences of means
##          dias tipoJugo
##         1.208    1.208
## replic.    12       12

Efectos

model.tables(x = m2,
             type = "effects",
             se = TRUE)

## Tables of effects
## 
##  dias 
## dias
##      0      3      7 
##  3.356 -0.678 -2.678 
## 
##  tipoJugo 
## tipoJugo
##       A       B       C 
##  1.0556  0.2056 -1.2611 
## 
## Standard errors of effects
##           dias tipoJugo
##         0.8542   0.8542
## replic.     12       12

Promedios e IC

ggeffect(model = m2)

## $dias
## # Predicted values of vitamC
## 
## dias | Predicted |         95% CI
## ---------------------------------
## 0    |     51.25 | [49.51, 52.99]
## 3    |     47.22 | [45.47, 48.96]
## 7    |     45.22 | [43.47, 46.96]
## 
## $tipoJugo
## # Predicted values of vitamC
## 
## tipoJugo | Predicted |         95% CI
## -------------------------------------
## A        |     48.95 | [47.21, 50.69]
## B        |     48.10 | [46.36, 49.84]
## C        |     46.63 | [44.89, 48.38]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "m2"

Comparaciones múltiples

Tukey

TukeyHSD(m2) %>% tidy()

## # A tibble: 6 × 7
##   term     contrast null.value estimate conf.low conf.high adj.p.value
##   <chr>    <chr>         <dbl>    <dbl>    <dbl>     <dbl>       <dbl>
## 1 dias     3-0               0   -4.03     -7.01    -1.06    0.00606  
## 2 dias     7-0               0   -6.03     -9.01    -3.06    0.0000634
## 3 dias     7-3               0   -2        -4.97     0.973   0.238    
## 4 tipoJugo B-A               0   -0.850    -3.82     2.12    0.763    
## 5 tipoJugo C-A               0   -2.32     -5.29     0.657   0.151    
## 6 tipoJugo C-B               0   -1.47     -4.44     1.51    0.454

Tamaño del efecto

eta_squared(m2)

## # Effect Size for ANOVA (Type I)
## 
## Parameter | Eta2 (partial) |       90% CI
## -----------------------------------------
## dias      |           0.46 | [0.22, 0.61]
## tipoJugo  |           0.11 | [0.00, 0.27]

Regresión Lineal Simple

Datos

cerezos_negros <- trees
cerezos_negros

##    Girth Height Volume
## 1    8.3     70   10.3
## 2    8.6     65   10.3
## 3    8.8     63   10.2
## 4   10.5     72   16.4
## 5   10.7     81   18.8
## 6   10.8     83   19.7
## 7   11.0     66   15.6
## 8   11.0     75   18.2
## 9   11.1     80   22.6
## 10  11.2     75   19.9
## 11  11.3     79   24.2
## 12  11.4     76   21.0
## 13  11.4     76   21.4
## 14  11.7     69   21.3
## 15  12.0     75   19.1
## 16  12.9     74   22.2
## 17  12.9     85   33.8
## 18  13.3     86   27.4
## 19  13.7     71   25.7
## 20  13.8     64   24.9
## 21  14.0     78   34.5
## 22  14.2     80   31.7
## 23  14.5     74   36.3
## 24  16.0     72   38.3
## 25  16.3     77   42.6
## 26  17.3     81   55.4
## 27  17.5     82   55.7
## 28  17.9     80   58.3
## 29  18.0     80   51.5
## 30  18.0     80   51.0
## 31  20.6     87   77.0

Objetivo y modelo

Objetivo: predecir el volumen de la madera en función del diámetro
Modelo: \(y_i = \beta_0 + \beta_1 \times X + \epsilon_i\), donde \(X: Diámetro\)

Gráfico

library(GGally)
ggpairs(cerezos_negros)

Nivel de significancia

En este caso usamos 5%.

Correlación

Hipótesis

\[H_0: \rho = 0 \\ H_0: \rho \neq 0 \]

Normalidad

Normalidad de diámetro:

shapiro.test(cerezos_negros$Girth)

## 
##  Shapiro-Wilk normality test
## 
## data:  cerezos_negros$Girth
## W = 0.94117, p-value = 0.08893

Normalidad de volumen:

shapiro.test(cerezos_negros$Volume)

## 
##  Shapiro-Wilk normality test
## 
## data:  cerezos_negros$Volume
## W = 0.88757, p-value = 0.003579

Spearman

cor.test(x = cerezos_negros$Girth, y = cerezos_negros$Volume, method = "spearman",
         conf.level = 0.95)

## 
##  Spearman's rank correlation rho
## 
## data:  cerezos_negros$Girth and cerezos_negros$Volume
## S = 224.61, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9547151

Modelo

modelo_regresion <- lm(Volume ~ Girth, data = cerezos_negros)
summary(modelo_regresion)

## 
## Call:
## lm(formula = Volume ~ Girth, data = cerezos_negros)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.065 -3.107  0.152  3.495  9.587 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -36.9435     3.3651  -10.98 7.62e-12 ***
## Girth         5.0659     0.2474   20.48  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.252 on 29 degrees of freedom
## Multiple R-squared:  0.9353, Adjusted R-squared:  0.9331 
## F-statistic: 419.4 on 1 and 29 DF,  p-value: < 2.2e-16

Residuales

residuales <- residuals(modelo_regresion)

Normalidad

shapiro.test(residuales)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.97889, p-value = 0.7811

Homocedasticidad

library(lmtest)
bptest(modelo_regresion)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_regresion
## BP = 5.6197, df = 1, p-value = 0.01776

Predicción

Estimación puntual

predict(object = modelo_regresion,
        newdata = data.frame(Girth = 20.1))

##        1 
## 64.88025

Estimación por intervalo

predict(object = modelo_regresion,
        newdata = data.frame(Girth = 20.1), interval = "confidence")

##        fit      lwr     upr
## 1 64.88025 61.07811 68.6824

Análisis estadístico de experimentos

Edimer David Jaramillo

2021-10-27

Experimento con 1 factor

Datos

Resumen descriptivo básico

Boxplot

Modelo estadístico

Hipótesis

Nivel de significancia

Anova

Análisis de residuales

Normalidad

Homocedasticidad

Estimaciones

Promedios

Efectos

Promedios e IC

Comparaciones múltiples

Sin ajuste

Con ajuste

Tukey

Reporte de resultados

Gráficos

Tamaño del efecto

DBCA

Datos

Estadísticos descriptivos

Por Tratamiento

Por Bloque

Gráfico

Modelo

Hipótesis

Principal

Verificación

Nivel de significancia

Anova

Análisis de residuales

Normalidad

Homocedasticidad

Estimaciones

Promedios

Efectos

Promedios e IC

Comparaciones múltiples

Tukey

Dunnet

Tamaño del efecto

Consecuencias de no bloquear

Modelo

Tukey

Diseño factorial 2 vías

Datos

Estadísticos descriptivos

Tipo de jugo

Días

Días y Tipo de jugo

Gráfico 1

Gráfico 2

Modelo

Hipótesis

Tipo de jugo

Días

Interacción

Nivel de significancia

Anova

Modelo 1

Modelo 2

Análisis de residuales

Normalidad

Homocedasticidad

Estimaciones

Promedios

Efectos

Promedios e IC

Comparaciones múltiples

Tukey

Tamaño del efecto

Regresión Lineal Simple

Datos