Last updated: 2021-06-12

Checks: 7 0

Knit directory: 2021/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20210412)

The command set.seed(20210412) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 520708a

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 520708a. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  Curso_Bioestadistica_MTripp_cuatriII.docx
    Untracked:  Curso_Bioestadistica_MTripp_cuatriII.pdf
    Untracked:  Diapositivas/
    Untracked:  README.html
    Untracked:  Resources/
    Untracked:  analysis/images/
    Untracked:  code/tarea_macrograd.R
    Untracked:  data/CS_subset.csv
    Untracked:  data/Consumo_oxigeno_wide.csv
    Untracked:  data/Darwin_esp.csv
    Untracked:  data/Data_enzimas_Experimento1.txt
    Untracked:  data/Data_enzimas_Experimento2.txt
    Untracked:  data/Data_enzimas_Experimento3.txt
    Untracked:  data/Data_enzimas_Experimento4.txt
    Untracked:  data/DownloadFestival(No Outlier).dat
    Untracked:  data/Festival.csv
    Untracked:  data/LungCapData.txt
    Untracked:  data/LungCapDataEsp.csv
    Untracked:  data/RExam.dat
    Untracked:  data/Rexamendat.csv
    Untracked:  data/Tabla1_Muestreo.txt
    Untracked:  data/Transcriptome_Anotacion.csv
    Untracked:  data/Transcriptome_DGE.csv
    Untracked:  data/datasets_Pokemon.csv
    Untracked:  data/datasets_Pokemon.xls
    Untracked:  data/exp_macrogard_growth.tab
    Untracked:  data/exp_macrogard_rna-dna.tab
    Untracked:  data/macrogard_crecimiento.csv
    Untracked:  data/penguins_size.csv
    Untracked:  data/pokemon_extended.csv
    Untracked:  output/Plot_all_penguins.pdf
    Untracked:  output/Plot_all_penguins.tiff
    Untracked:  output/graficos/

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/Clase7_tstudent.Rmd) and HTML (docs/Clase7_tstudent.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	520708a	Miguel Tripp	2021-06-12	Publish the initial files for myproject
html	99c3644	Miguel Tripp	2021-06-05	Build site.
Rmd	624db7c	Miguel Tripp	2021-06-05	Publish the initial files for myproject
html	fbd7857	Miguel Tripp	2021-06-05	Build site.
Rmd	5028622	Miguel Tripp	2021-06-05	Publish the initial files for myproject
html	87a646f	Miguel Tripp	2021-06-02	Build site.
Rmd	5e69a71	Miguel Tripp	2021-06-02	Publish the initial files for myproject

library(webr)
library(patchwork)
library(tidyverse)

-- Attaching packages --------------------------------------- tidyverse 1.3.1 --

v ggplot2 3.3.3     v purrr   0.3.4
v tibble  3.0.4     v dplyr   1.0.5
v tidyr   1.1.3     v stringr 1.4.0
v readr   1.4.0     v forcats 0.5.1

-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

1 Distribución T

La distribución T-student se asemeja en gran medida a la distribución normal. Tiene como parámetros la media, la varianza y además incorpora a través de los grados de libertad una modificación que permite flexibilizar las colas en función del tamaño que tenga la muestra. A medida que se reduce el tamaño muestral, la probabilidad acumulada en las colas aumenta, siendo así menos estricta de lo cabría esperar en una distribución normal. Una distribución T-student con 30 o más grados de libertad es prácticamente igual a una distribución normal.

null device 
          1

2 Comparación de dos medias con t de Student

La prueba t estima la diferencia entre la media de dos grupos usando la proporción de la diferencía entre la media de los grupos y el error estandar de ambos grupos

2.1 Supuestos de la prueba

Para esta prueba es necesario considerar el nivel de medición de las variables, para lo cual se requiere de una variable cuantitativa y una variable nominal (dos grupos).

Las condiciones para calcular intervalos de confianza o aplicar un test basados en la distribución T-student son:

Independencia: Las observaciones tienen que ser independientes unas de las otras. Para ello el muestreo tiene que ser aleatoria y el tamaño de la muestra inferior al 10% de la población.
Normalidad: Las poblaciones que se comparan tienn que distribuirse de forma normal. En caso de cierta asimetría, los T-test son considerablemente robustos cuando el tamaño de muestra es mayor o igual a 30.
Homogeneidad de varianzas: La varianza de ambas poblaciones comparadas debe ser igual. En caso de no cumplirse con esta condición se puede usar una prueba Welch Two Sample t-test. Esta corrección se incorpora a traves de los grados de libertad permitiendo compensar la diferencia de varianzas (con el inconveniente de que se pierde precisión).

2.2 Establecer las hipótesis

La hopótesis nula (\(H_0\)) considera que no hay diferencia o cambio. Suele contener en su definición el simbolo \(=\) de manera que al comparar dos medias independientes, suponemos que \(\mu_1 = \mu_2\). La hipotesis alternativa (\(H_A\)) considera que el valor real de las medias de ambas poblaciones son diferentes, de manera que \(\mu_1 \ne \mu_2\)

2.3 Determinar el tipo de test; una o dos colas

Los test de hipótesis pueden ser de una cola o de dos colas. Si la hipótesis alternativa emplea > o < se trata de un test de una cola, en el que solo se analizan desviaciones en un sentido. Si la hipótesis alternativa es del tipo “diferente de” se trata de un test de dos colas, en el que se analizan posibles desviaciones en las dos direcciones. Solo se emplean test de una cola cuando se sabe con seguridad que las desviaciones de interés son en un sentido y únicamente si se ha determinado antes de observar la muestra, no a posteriori.

2.4 Formula

\(t = \frac{\tilde{x}_1 - \tilde{x}_2} {\sqrt{ \frac{S^2_{1}}{n_{1}} + \frac{S^2_{2}}{n_{2}} }}\)

Donde t es es valor t, \(\tilde{x}_1\) y \(\tilde{x}_2\) son las mediates de los grupos contrastados, \(S_1\) y \(S_2\) son las desviaciones estandar de cada grupo, y \(n_1\) y \(n_2\) es el número de observaciones en cada grupo.

Un valor de t grande indicaria que la diferencia entre la media de los grupos es mayor que la desviación estandar de ambos grupos, indicando una diferencia entre ambos grupos.

Posteriormente comparamos el valor t calculado contra valores de probabilidad.

El nivel de significancia \(\alpha\) determina la probabilidad de error que se quiere asumir a la hora de rechazar la hipótesis nula. Se emplea como punto de referencia para determinar si el valor p obtenido en el test de hipótesis es suficientemente bajo como para considerar significativas las diferencias observadas y por lo tanto rechazar \(H_0\).

3 Ejemplo (manual):

Se quiere estudiar el efecto de un aditivo en el consumo de alimento en una dieta de peces, para ello, se mide la ingesta de comida (gr) en dos grupos de peces. En vista de los resultados ¿se puede conisderar que el aditivo funciona en la tasa de ingesta?

Resultados del experimento
Grupo	Media	S	n
control	24.889	8.007	22
tratamiento	38.23	15.398	22

1. Planteamos las hipótesis

\(H_0:\) No hay diferencias entre las medias; media control - media tratamiento \(= 0\)

\(H_a:\) Hay diferencias entre las medias; media control - media tratamiento \(\ne 0\)

2. Páramtetro estimado

media control - media tratamiento = 24.889 - 38.23 = -13.341

4. Modelo de la prueba

Como queremos evaluar si hay diferencias por las dietas (que puede ser positiva o negativa) usamos una prueba de dos colas

5. Nivel de significancia

\(\alpha = 0.05\)

6. Estimación del valor p

Parametro estimado = -13.341

Grados de libertad = 21

SE(media control - media tratamiento) =

\(SEM = \sqrt{ \frac{S^2_{Control}}{n_{control}} + \frac{S^2_{tratamiento}}{n_{tratamiento}} }\)

Si lo reemplazamos:

SEM = sqrt((8.007^2 / 22) + (15.398^2 / 22))

[1] 3.700187

Estimación del valor T

T = 38.23 - 24.889 / 3.7001871

[1] -3.605493

Por último, estimamos el valor de probabilidad

2*pt(-abs(T), df=42)

[1] 0.0008203744

Para aplicar la prueba de t de Student en R, es posible utilizar la función t.test() la cual esta incluida en R base y utiliza la siguiente syntaxis:

t.test(x, y , alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
t.test(variable1 ~ variable2, data = datos, var.equal = TRUE)

en donde var.equal= nosotros indicamos si ambos grupos muestreales tienen varianzas homogeneas o no.

t.test(ctrl, treat, var.equal = TRUE)


    Two Sample t-test

data:  ctrl and treat
t = -3.6055, df = 42, p-value = 0.0008203
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -20.807928  -5.873726
sample estimates:
mean of x mean of y 
 24.88915  38.22998

El resultadod e esta prueba indica:

Idicación de que es lo que se esta comparando (data).
El valor t. Es frecuente obtener valores negativos, lo cual esta bien, ya que solo nos interesa el valor absoluto de la diferencía, o la idstancia de 0, sin importat la dirección.
Los grados de libertad, los cuales estan relacionados con el tamaño de la muestra e indican cuantos datos “lilbres” estan disponibles para hacer las comparaciones.
El valor p el cual describe la probabilidad de ver un valor t tan extremo o mas como el obtenido por azar.
La hipotesis alternativa (\(H_a\)). En este caso, la \(H_a\) es que la diferencía entre la media de las muestras no es 0.
Intervalo de confianza del 95%. Esto significa que de repetir el muestreo multiples veces, el 95% de las veces la media estaraá contenida en este intervalo.
La media calculada para cada grupo

4 Prueba T para una sola muestra

Esta prueba se utiliza cuando se comara la media de una muestra contra una media (\(\mu\)) conocida o hipotética

Generalmente, este valor teórico de la media puee venir de un experimento previo o de un valor abitrario.

Para realizar la prueba t de una muestra en R se puede utilizar la función t.test() con la siguiente estructura:

t.test(x, y = NULL,
       alternative = c("two.sided", "less", "greater"),
       mu = 0, paired = FALSE, var.equal = FALSE,
       conf.level = 0.95, ...)

en donde:

x: vector numérico con los datos alternative: tipo de hipotesis alterna. Los valores disponibles son two.sided cuando la hipótesis alternativa es \(x \ne \mu\), less para el caso \(x < \mu\) y greater para el caso \(x > \mu\).

4.1 Ejercicio:

Para verificar si el proceso de llenado de bolsa de un alimento seco para cultivo de 500 g esta operando correctamente, se toman aleatoriamente muestras de tamaño de diez bolsas cada 4 horas Una de estas muestars esta compuesta por las siguientes observaciones:

510, 492, 494, 498, 492, 496, 502, 491, 507, 496

Entonces, vamos a responder la pregunta: a un nivel de signifcancia del 5% ¿esta el proceso de llenado de la bolsa llevandose a cabo correctamente?

El primer paso consiste en explorar si la muestra viene de una distribución normal.

bolsa <- c(510, 492, 494, 498, 492, 496, 502, 491, 507, 496)

# Prueba de normalidad shapiro
shapiro.test(bolsa)


    Shapiro-Wilk normality test

data:  bolsa
W = 0.88468, p-value = 0.1476

Como el valor P de la prueba de Shapiro-Wilk es mayor a 0.05, se puede asumir que la muestra de bolsas proviene de una distribución normal.

Plantear las hipótesis

\(H_0: \mu = 500 gr\)

\(H_A: \mu \ne 500gr\)

Realizar la prueba

t.test(bolsa, alternative = "two.sided", mu = 500)


    One Sample t-test

data:  bolsa
t = -1.0629, df = 9, p-value = 0.3155
alternative hypothesis: true mean is not equal to 500
95 percent confidence interval:
 493.1176 502.4824
sample estimates:
mean of x 
    497.8

plot(t.test(bolsa, alternative = "two.sided", mu = 500))

Como el vaor P es mayor al 5% entonces concluimos que no hay evidencia para afirmar que el proceso de llenao no se esta realizando de manera incorrecta.

5 Comparación de dos medias independientes

5.1 Ejemplo 1: ¿Hay diferencia entre la longitud del culmen entre las poblaciones de pingünos?

De la base de datos de pingüinos, queremos ver si hay diferencias en la longitud del culmen entre entre las poblaciones de Dream y Biscoe de la especie Adelie

Para esta sesión, necesitaremos tres paquetes:

Tidyverse: Para usar dplyr y ggplot2
car: contiene algunas funciones estadisticas
rstatixs: Contiene funciones estadisticas compatibles con el ambiente de tidyverse (pipe-friendly). Puedes conocer mas sobre este paquete aqui

library(tidyverse)
library(car)
library(rstatix)

penguins <- read_csv("data/penguins_size.csv")


# Filtrar especies
pen_Adelie <- penguins %>%
  filter(species == "Adelie", island %in% c("Dream", "Biscoe"))

El primer paso es realizar una exploración de los datos y verificar los supuestos de la prueba; la prueba t requiere que los datos tengan una distribución normal y con varianzas iguales.

Grafico de histogramas para visualizar la distribución de los datos entre ambas islas

ggplot(pen_Adelie, aes(x = culmen_length_mm, fill = island, col = island))+
  geom_histogram(alpha = 0.3)+
  facet_wrap(~island)+
  labs(y = "frecuencia",
       x = "longitud del culmen (mm)",
       title="Grafico de densidad", 
       fill = "isla", 
       col = "isla")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Tambien es posible realizar un gráfico qq (quantil-quantil) para evaluar visualmente si la distribución de los datos de ambas islas se asemejan a una distribución normal

ggplot(pen_Adelie, aes(sample = culmen_length_mm, col = island))+
  stat_qq()+
  stat_qq_line()+
  facet_grid(.~ island)

De acuerdo a los gráficos, es evidente que ambas poblaciones tienen distribuciones muy parecida a la normal. Para validar estos resutlados, realizaremos una prueba de Shapiro_wilks para evaluar la significancia de la normalidad.

El método de Shapiro-Wilks es ampliamente recomendado y tiene mayor poder en comparación de K-S. Para conocer mas sobre pruebas de normalidad visita aqui

pin_Ade_shap <- pen_Adelie %>%
  group_by(island) %>%
  shapiro_test(culmen_length_mm) %>%
  print()

# A tibble: 2 x 4
  island variable         statistic     p
  <chr>  <chr>                <dbl> <dbl>
1 Biscoe culmen_length_mm     0.977 0.520
2 Dream  culmen_length_mm     0.984 0.640

La hipótesis nula de esta prueba es que la muestra viene de una población con distribución normal.

Tanto la gráfica de densidad como el qqplot sugieren que los datos cumplen con el supuesto de normalidad. Para validar esta aseveración es necesario realizar la prueba de Shapiro-Wilk.

El siguiente paso es evaluar el supuesto de homogeneidad de varianzas. Las gráficas de densidades nuevamente nos dan un indicativo visual de la disperción de los datos, pero es posible evaluar la significancia con una prueba de Levene implementada en el paquete car:

leveneTest(culmen_length_mm ~ island, data = pen_Adelie)

Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.0054 0.9415
      98

pen_Adelie %>% 
  levene_test(culmen_length_mm ~ island)

Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.

# A tibble: 1 x 4
    df1   df2 statistic     p
  <int> <int>     <dbl> <dbl>
1     1    98   0.00542 0.941

La prueba nos dice que las varianzas son homogeneas entre los grupos

El siguiente paso para la exploración de los datos es obtener el promedio de los valores de la longitud del culmen:

pen_Ade_sum <- pen_Adelie %>%
  group_by(island) %>%
  summarise(promedio = mean(culmen_length_mm),
            desviacion = sd(culmen_length_mm),
            N = n())
pen_Ade_sum

# A tibble: 2 x 4
  island promedio desviacion     N
  <chr>     <dbl>      <dbl> <int>
1 Biscoe     39.0       2.48    44
2 Dream      38.5       2.47    56

De manera alternativa, este procedimiento se puede realizar con ayuda del paquete rstatix

pen_Ade_sum <- pen_Adelie %>%
  group_by(island) %>%
  get_summary_stats(culmen_length_mm)
pen_Ade_sum

# A tibble: 2 x 14
  island variable       n   min   max median    q1    q3   iqr   mad  mean    sd
  <chr>  <chr>      <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Biscoe culmen_le~    44  34.5  45.6   38.7  37.7  40.7  3.02  2.74  39.0  2.48
2 Dream  culmen_le~    56  32.1  44.1   38.6  36.8  40.4  3.62  2.74  38.5  2.46
# ... with 2 more variables: se <dbl>, ci <dbl>

Posteriormente, hacemos un boxplot para visualizar los datos:

ggplot(pen_Adelie, aes(x = island, y = culmen_length_mm, fill = island))+
  geom_boxplot(outlier.shape = "")+
  geom_point(position = position_jitter(0.1), col = "grey50")+
  theme_light()+
  labs(x = "Isla",
       y = "longitud del culmen (mm)",
       title = "Longitud del culmen por isla en el pingüino Adelie",
       fill = "Isla", col = "Isla")

Finalmente, hacemos la prueba t de student. Esta puede realizarse siguiendo el comando base en R t.test()

t.test(culmen_length_mm ~ island, data = pen_Adelie, var.equal = TRUE )


    Two Sample t-test

data:  culmen_length_mm by island
t = 0.95016, df = 98, p-value = 0.3444
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.5151265  1.4615551
sample estimates:
mean in group Biscoe  mean in group Dream 
            38.97500             38.50179

En esta prueba, la hipótesis nula es que la media de las dos poblaciones es igual, lo que impricla que la hipótesis alternativa es que la diferencia de las medias no es igual a cero.

En la prueba de hiótesis, se desea aceptar o rechazar la hipotesis nula con ciertos intervalos de confianza. Dado que probamos la diferencias entre ambas medias, el intervalo de confianza especifica el intervalo de los valores de esta diferencia.

Alternativametne podemos usar la prueba implementada en rstatixs:

pin_Ade_ttest <- pen_Adelie %>%
  t_test(culmen_length_mm ~ island, var.equal = TRUE) %>%
  print()

# A tibble: 1 x 8
  .y.              group1 group2    n1    n2 statistic    df     p
* <chr>            <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
1 culmen_length_mm Biscoe Dream     44    56     0.950    98 0.344

Ahora con la información generada con la prueba, es posible agregar la anotacion a la gráfica

p.val = pin_Ade_ttest$p

ggplot(pen_Adelie, aes(x = island, y = culmen_length_mm, fill = island))+
  geom_boxplot(outlier.shape = "")+
  geom_point(position = position_jitter(0.1), col = "grey50")+
  theme_light()+
  labs(x = "Isla",
       y = "longitud del culmen (mm)",
       title = "Longitud del culmen por isla en el pingüino Adelie",
       fill = "Isla", col = "Isla",
       caption = paste0("t-student; p = ", p.val))

Tambien es posible visualizar la ditrtibución de probabilidad

plot(t.test(culmen_length_mm ~ island, data = pen_Adelie, var.equal = TRUE ))

5.2 Ejemplo 2: ¿ Hay diferencia en la longitud de la aleta entre pingünio Adelie y Gentoo?

pen_Aleta <- penguins %>%
  filter(species %in% c("Adelie", "Gentoo"))

Histograma de frecuencias

ggplot(pen_Aleta, aes(x = flipper_length_mm, fill = species, col = species))+
  geom_histogram(alpha = 0.3)+
  facet_wrap(~species)+
  labs(y = "Frecuencia",
       x = "longitud de la aleta (mm)",
       title="Grafico de densidad", 
       fill = "especie", 
       col = "especie")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 2 rows containing non-finite values (stat_bin).

qqplot

ggplot(pen_Aleta, aes(sample = flipper_length_mm, col = species))+
  stat_qq()+
  stat_qq_line()+
  facet_grid(.~ species)

Warning: Removed 2 rows containing non-finite values (stat_qq).

Warning: Removed 2 rows containing non-finite values (stat_qq_line).

Prueba de shapiro utilizando el paquete RSTATIX

pen_Aleta_shap <- pen_Aleta %>%
  group_by(species) %>%
  shapiro_test(flipper_length_mm) %>%
  print()

# A tibble: 2 x 4
  species variable          statistic       p
  <chr>   <chr>                 <dbl>   <dbl>
1 Adelie  flipper_length_mm     0.993 0.720  
2 Gentoo  flipper_length_mm     0.962 0.00162

Distribucion de las varianzas

leveneTest(flipper_length_mm ~ species, data = pen_Aleta)

Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.

Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   1   0.157 0.6922
      272

Boxplot

ggplot(pen_Aleta, aes(x = species, y = flipper_length_mm, fill = species))+
  geom_boxplot()+
  geom_point(position = position_jitter(0.1), col = "grey50")+
  theme_light()+
  labs(x = "Especie",
       y = "longitud de la aleta (mm)",
       title = "Longitud de la aleta por especie",
       fill = "Especie", col = "Especie")

Warning: Removed 2 rows containing non-finite values (stat_boxplot).

Warning: Removed 2 rows containing missing values (geom_point).

Hacemos la prueba T de student

t.test(flipper_length_mm ~ species, data = pen_Aleta, var.equal = TRUE)


    Two Sample t-test

data:  flipper_length_mm by species
t = -34.415, df = 272, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -28.79125 -25.67545
sample estimates:
mean in group Adelie mean in group Gentoo 
            189.9536             217.1870

t.student con rstatixs

pen_Aleta_ttest <- pen_Aleta %>%
  t_test(flipper_length_mm ~ species, var.equal = TRUE)
pen_Aleta_ttest

# A tibble: 1 x 8
  .y.               group1 group2    n1    n2 statistic    df         p
* <chr>             <chr>  <chr>  <int> <int>     <dbl> <dbl>     <dbl>
1 flipper_length_mm Adelie Gentoo   152   124     -34.4   272 4.21e-101

La prueba nos dice que las medias de ambas poblaciones son diferentes.

5.3 Ejemplo 3: ¿Hay diferencia en el peso de los pokemon de roca y los de planta?

Para resolver este ejercicio, vamos a hacer uso de la tabla datasets_pokemon.xlsy pokemon_extended.csv para relacionar el tipo y el peso de cada pokemon.

poke_ext <- read_csv("data/pokemon_extended.csv")


-- Column specification --------------------------------------------------------
cols(
  height_m = col_double(),
  Name = col_character(),
  percentage_male = col_double(),
  pokedex_number = col_double(),
  weight_kg = col_double()
)

pokemon <- readxl::read_xls("data/datasets_pokemon.xls")


poke_test <- pokemon %>% 
  inner_join(poke_ext, by = "Name") %>% 
  select(Type1, weight_kg) %>% 
  filter(Type1 %in% c("Grass", "Rock"))

Como primer paso, visualizamos los valores de peso de cada grupo con un histograma

poke_hito <- ggplot(poke_test, aes(x = weight_kg, fill = Type1))+
  geom_histogram(alpha = 0.4)+
  facet_wrap(~Type1, nrow = 2)+
  theme(legend.position = "none")

y psoteriorment construimos qqplots para cada grupo

poke_qq <- ggplot(poke_test, aes(sample = weight_kg, color = Type1))+
  geom_qq()+
  geom_qq_line()+
  facet_wrap(~Type1, nrow = 2)+
  theme(legend.position = "none")


poke_hito + poke_qq

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 4 rows containing non-finite values (stat_bin).

Warning: Removed 4 rows containing non-finite values (stat_qq).

Warning: Removed 4 rows containing non-finite values (stat_qq_line).

Los histogramas muestran que los valores parecen tener una fuerte asimetria y curtosis negativa para el caso del grupo roca.

Prueba de Shapiro-Wilk

poke_test %>% 
  group_by(Type1) %>% 
  shapiro_test(weight_kg)

# A tibble: 2 x 4
  Type1 variable  statistic            p
  <chr> <chr>         <dbl>        <dbl>
1 Grass weight_kg     0.711 0.0000000429
2 Rock  weight_kg     0.758 0.0000871

La prueba de Shapiro-Wilk nos confirma que ambos grupos no vienen de una distribución normal.En este caso, es posible hacer una transformación de los datos, para lo cual vamos a generar una columna columna con el logaritmo del peso

poke_test$log_weight <- log(poke_test$weight_kg)
poke_test

# A tibble: 72 x 3
   Type1 weight_kg log_weight
   <chr>     <dbl>      <dbl>
 1 Grass       6.9       1.93
 2 Grass      13         2.56
 3 Grass     100         4.61
 4 Grass       5.4       1.69
 5 Grass       8.6       2.15
 6 Grass      18.6       2.92
 7 Grass       4         1.39
 8 Grass       6.4       1.86
 9 Grass      15.5       2.74
10 Rock       NA        NA   
# ... with 62 more rows

Y volvemos a generar los histogramas de frecuencia y qqplots para ser el cambio

poke_hito_log <- ggplot(poke_test, aes(x = log_weight, fill = Type1))+
  geom_histogram(alpha = 0.4)+
  facet_wrap(~Type1, nrow = 2)+
  theme(legend.position = "none")

poke_qq_log <- ggplot(poke_test, aes(sample = log_weight, color = Type1))+
  geom_qq()+
  geom_qq_line()+
  facet_wrap(~Type1, nrow = 2)+
  theme(legend.position = "none")


poke_hito_log + poke_qq_log

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 4 rows containing non-finite values (stat_bin).

Warning: Removed 4 rows containing non-finite values (stat_qq).

Warning: Removed 4 rows containing non-finite values (stat_qq_line).

y repetimos la prueba de normalidad

poke_test %>% 
  group_by(Type1) %>% 
  shapiro_test(log_weight)

# A tibble: 2 x 4
  Type1 variable   statistic     p
  <chr> <chr>          <dbl> <dbl>
1 Grass log_weight     0.972 0.349
2 Rock  log_weight     0.930 0.110

Comprobando que la transformación nos ajusta los datos a una distribución normal.

Posteriormente, realizamos una prueba de Levene para corroborar que los datos tengan homogeneidad de varianzas

poke_test %>% 
  levene_test(log_weight ~ Type1)

Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.

# A tibble: 1 x 4
    df1   df2 statistic     p
  <int> <int>     <dbl> <dbl>
1     1    66     0.239 0.627

Posteriormente realizamos la prueba t usando los valores transformados

t.test(log_weight ~ Type1, data = poke_test)


    Welch Two Sample t-test

data:  log_weight by Type1
t = -3.3259, df = 49.745, p-value = 0.001662
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.8627433 -0.4598902
sample estimates:
mean in group Grass  mean in group Rock 
           2.624765            3.786082

plot(t.test(log_weight ~ Type1, data = poke_test))

Finalmente, realizamos un boxplot para mostrar las diferencias entre ambos grupos

ggplot(poke_test, aes(x = Type1, y = weight_kg, fill = Type1))+
  geom_boxplot()

Warning: Removed 4 rows containing non-finite values (stat_boxplot).

ggplot(poke_test, aes(x = Type1, y = log_weight, fill = Type1))+
  geom_boxplot()

Warning: Removed 4 rows containing non-finite values (stat_boxplot).

5.4 Ejemplo 4: ¿Hay diferencia entre la actividad enzimática de la citrato sintasa a 18°C y 32°C?; Prueba t con varianzas heterogeneas

Como ya hemos visto, una de las condiciones para aplicar una prueba t de Student es la igualdad de varianzas (homocedasticidad) de ambas poblaciones. En caso de no cumplirse esta condición se puede emplear un Welch Two Sample t-test. Esta corrección se incorpora a través de los grados de libertad permitiendo compensar la diferencia de varianzas. El número de grados de libertad de un Welch Two Sample t-test viene dado por la siguiente función:

\(f = \frac{(\frac{S^2_1}{n_1} + \frac{S^2_2}{n_2})^2} {\frac{1}{n_1 + 1}(\frac{S^2_1}{n_1})^2 + \frac{1}{n_2 + 1}(\frac{S^2_2}{n_2})^2} - 2\)

Para este ejercicio, utilizaremos los datos de actividad enzimática de la citrato sintasa (CS) entre una muestra colectada tras una exposición a 18°C y una exposición a 32°C. Estos datos se pueden descargar de la siguiente liga:

cs_url <- "https://raw.githubusercontent.com/trippv/Miguel_Tripp/master/CS_subset.csv"

cs_activity <- read_csv(cs_url)

Para este ejercicio vamos a filtrar solamente los datos del experimento número 4 y de las temperaturas 18 y 32. Posteriomente vamos a convertir estos valores a factores

cs_activity <- cs_activity %>% 
  filter(ExpNum== 4,
         TreatTemp %in% c(18,32)) %>% 
  mutate(TreatTemp = factor(TreatTemp))

Utilizando rstatix, realizamos las pruebas de normalidad y homogeneidad de varianzas

# Prueba de normalidad
cs_activity %>% 
  group_by(TreatTemp) %>% 
  shapiro_test(ActivityCS)

# A tibble: 2 x 4
  TreatTemp variable   statistic      p
  <fct>     <chr>          <dbl>  <dbl>
1 18        ActivityCS     0.882 0.0930
2 32        ActivityCS     0.955 0.704

# Prueba de homogeneidad de varianzas
cs_activity %>% 
  levene_test(ActivityCS ~ TreatTemp)

# A tibble: 1 x 4
    df1   df2 statistic       p
  <int> <int>     <dbl>   <dbl>
1     1    22      10.5 0.00377

# Boxplot 

ggplot(data = cs_activity, aes(x = TreatTemp, y = ActivityCS, fill = TreatTemp))+
  geom_boxplot()

Finalmente realizamos la prueba de Welch incluyendo el parametro var.equal = FALSE

cs_activity %>% 
  t_test(ActivityCS ~ TreatTemp, var.equal = FALSE)

# A tibble: 1 x 8
  .y.        group1 group2    n1    n2 statistic    df      p
* <chr>      <chr>  <chr>  <int> <int>     <dbl> <dbl>  <dbl>
1 ActivityCS 18     32        12    12      1.99  16.3 0.0639

# alternativamente
t.test(ActivityCS ~ TreatTemp, data = cs_activity, var.equal = FALSE)


    Welch Two Sample t-test

data:  ActivityCS by TreatTemp
t = 1.9884, df = 16.275, p-value = 0.06385
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.007403253  0.236401539
sample estimates:
mean in group 18 mean in group 32 
       0.6242346        0.5097355

plot(t.test(ActivityCS ~ TreatTemp, data = cs_activity, var.equal = FALSE))

Finalmente concluimos que no hay diferencia entre la actividad enzimatica tras una exposición a 32°C

6 Comparación de dos medias dependientes

6.1 Ejercicio: ¿Hay diferencia en el consumo de oxígeno (\(MO_2\)) tras la exposición a la temperatura?

Para esto utilizaremos datos de tasa metabolica de 10 juveniles de abulón azul a 18°C y posteriormente medido tras una exposición aguda a 32°C, es decir cada invdividuo aparece en ambos grupos

abalone_mo <- "https://raw.githubusercontent.com/trippv/Miguel_Tripp/master/abalone_MO2.csv"

mo2 <- read_csv(abalone_mo)


-- Column specification --------------------------------------------------------
cols(
  Individuo = col_character(),
  Temperatura = col_double(),
  MO2 = col_double()
)

mo2$Temperatura <- as.factor(mo2$Temperatura)

Exploración de los datos con geom_histogram

ggplot(mo2, aes(x = MO2,fill=Temperatura))+
  geom_histogram()+
  labs(x = "consumo de oxígeno (MO2)")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Verificar la normalidad con gráfico de qq seguido con una prueba de Shapiro

ggplot(mo2, aes(sample= MO2, col = Temperatura)) +
  stat_qq()+
  stat_qq_line()+
  facet_grid(.~ Temperatura)

mo2_shapiro <- mo2 %>%
  group_by(Temperatura) %>%
  shapiro_test(MO2) %>% 
  print()

# A tibble: 2 x 4
  Temperatura variable statistic      p
  <fct>       <chr>        <dbl>  <dbl>
1 18          MO2          0.829 0.0329
2 32          MO2          0.879 0.128

Ambas gráficas y la prueba indica que uno de los grupos no cumple con el supuesto de normalidad. Dado que el tamaño de la muestra es pequeño, hay que evaluar cual es el mejor procedimiento.

Una opción es analizar cada una de las muestras y evaluar si alguno puede ser considerado como outlier. Otra opcion es realizar una transformación de los datos.

A continuación, se tranformaran los datos de MO2 a logaritmo

mo2 <- mo2 %>% 
  mutate(mo2_log = log(MO2))

Verificar la normalidad de los datos transformados con gráfico de qq seguido con una prueba de Shapiro

ggplot(mo2, aes(sample= mo2_log, col = Temperatura)) +
  stat_qq()+
  stat_qq_line()+
  facet_grid(.~ Temperatura)+
  labs(title = "QQplot de datos transformados")

Prueba de Shapiro para los datos transformados

mo2_log_shapiro <- mo2 %>%
  group_by(Temperatura) %>%
  shapiro_test(mo2_log) %>% 
  print()

# A tibble: 2 x 4
  Temperatura variable statistic     p
  <fct>       <chr>        <dbl> <dbl>
1 18          mo2_log      0.952 0.696
2 32          mo2_log      0.929 0.441

Prueba de levene para varianzas

mo2_log_leven <- mo2 %>%
  levene_test(mo2_log ~ Temperatura) %>% 
  print()

# A tibble: 1 x 4
    df1   df2 statistic      p
  <int> <int>     <dbl>  <dbl>
1     1    18      4.42 0.0499

Prueba t de student pareada

mo2_log_ttest <- mo2 %>%
  t_test(mo2_log ~ Temperatura, paired = TRUE, var.equal = TRUE)

Y finalmente graficamos los datos

ggplot(mo2, aes(x = Temperatura, y = MO2, fill = Temperatura))+
  geom_boxplot(outlier.shape = "")+
  geom_point(col = "grey35")+
  geom_line(aes(group = Individuo), col = "grey35")+
  labs(title = "Cambio en el consumo de oxígeno tras un incremento en la temperatura")

7 Reportar los resultados de una prueba t con `apa`

El paquete apa genera un reporte con los resultados de la prueba t de Student de acuerdo con las normas APA. Para esto se requiere un objeto t.test.

#install.packages("apa")
library(apa)


Attaching package: 'apa'

The following objects are masked from 'package:rstatix':

    cohens_d, t_test

# Ejemplo
prueba_t <- t.test(mo2_log ~ Temperatura, data = mo2, var.equal = TRUE, paired = FALSE)

t_apa(prueba_t)

t(18) = -8.42, p < .001, d = -3.77

sessionInfo()

R version 4.0.5 (2021-03-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] apa_0.3.3       rstatix_0.7.0   car_3.0-10      carData_3.0-4  
 [5] forcats_0.5.1   stringr_1.4.0   dplyr_1.0.5     purrr_0.3.4    
 [9] readr_1.4.0     tidyr_1.1.3     tibble_3.0.4    ggplot2_3.3.3  
[13] tidyverse_1.3.1 patchwork_1.1.1 webr_0.1.6      workflowr_1.6.2

loaded via a namespace (and not attached):
  [1] rrtable_0.2.1      colorspace_2.0-0   ggsignif_0.6.0    
  [4] ellipsis_0.3.1     rio_0.5.16         sjlabelled_1.1.8  
  [7] rprojroot_2.0.2    flextable_0.6.6    base64enc_0.1-3   
 [10] fs_1.5.0           rstudioapi_0.13    ggpubr_0.4.0      
 [13] farver_2.0.3       DT_0.17            fansi_0.4.2       
 [16] lubridate_1.7.10   xml2_1.3.2         mnormt_2.0.2      
 [19] knitr_1.30         sjmisc_2.8.7       polyclip_1.10-0   
 [22] jsonlite_1.7.2     broom_0.7.6        dbplyr_2.1.1      
 [25] ggforce_0.3.2      shiny_1.5.0        compiler_4.0.5    
 [28] httr_1.4.2         MBESS_4.8.0        backports_1.2.1   
 [31] assertthat_0.2.1   fastmap_1.0.1      cli_2.5.0         
 [34] later_1.1.0.1      tweenr_1.0.1       moonBook_0.2.4    
 [37] htmltools_0.5.1.1  tools_4.0.5        gtable_0.3.0      
 [40] glue_1.4.2         Rcpp_1.0.5         cellranger_1.1.0  
 [43] vctrs_0.3.8        nlme_3.1-152       psych_2.0.12      
 [46] lmtest_0.9-38      insight_0.11.1     xfun_0.20         
 [49] ps_1.5.0           openxlsx_4.2.3     rvest_1.0.0       
 [52] mime_0.9           miniUI_0.1.1.1     lifecycle_1.0.0   
 [55] devEMF_4.0-2       MASS_7.3-53        zoo_1.8-8         
 [58] scales_1.1.1       hms_1.0.0          promises_1.1.1    
 [61] parallel_4.0.5     RColorBrewer_1.1-2 yaml_2.2.1        
 [64] curl_4.3           gdtools_0.2.3      stringi_1.5.3     
 [67] zip_2.1.1          rlang_0.4.11       pkgconfig_2.0.3   
 [70] systemfonts_0.3.2  evaluate_0.14      lattice_0.20-41   
 [73] labeling_0.4.2     htmlwidgets_1.5.3  rvg_0.2.5         
 [76] tidyselect_1.1.1   editData_0.1.8     magrittr_2.0.1    
 [79] R6_2.5.0           generics_0.1.0     DBI_1.1.0         
 [82] pillar_1.6.0       haven_2.3.1        whisker_0.4       
 [85] foreign_0.8-81     withr_2.4.2        abind_1.4-5       
 [88] ztable_0.2.2       modelr_0.1.8       crayon_1.4.1      
 [91] shinyWidgets_0.5.4 uuid_0.1-4         utf8_1.2.1        
 [94] tmvnsim_1.0-2      rmarkdown_2.6      officer_0.3.18    
 [97] grid_4.0.5         readxl_1.3.1       data.table_1.13.6 
[100] git2r_0.27.1       vcd_1.4-8          reprex_2.0.0      
[103] digest_0.6.27      xtable_1.8-4       httpuv_1.5.4      
[106] munsell_0.5.0

Comparación de dos medias: Prueba de de Student

Miguel Tripp

2021-05-29

1 Distribución T

2 Comparación de dos medias con t de Student

2.1 Supuestos de la prueba

2.2 Establecer las hipótesis

2.3 Determinar el tipo de test; una o dos colas

2.4 Formula

3 Ejemplo (manual):

4 Prueba T para una sola muestra

4.1 Ejercicio:

5 Comparación de dos medias independientes

5.1 Ejemplo 1: ¿Hay diferencia entre la longitud del culmen entre las poblaciones de pingünos?

5.2 Ejemplo 2: ¿ Hay diferencia en la longitud de la aleta entre pingünio Adelie y Gentoo?

5.3 Ejemplo 3: ¿Hay diferencia en el peso de los pokemon de roca y los de planta?

5.4 Ejemplo 4: ¿Hay diferencia entre la actividad enzimática de la citrato sintasa a 18°C y 32°C?; Prueba t con varianzas heterogeneas

6 Comparación de dos medias dependientes

6.1 Ejercicio: ¿Hay diferencia en el consumo de oxígeno (\(MO_2\)) tras la exposición a la temperatura?

7 Reportar los resultados de una prueba t con `apa`

Comparación de dos medias: Prueba de de Student

Miguel Tripp

2021-05-29

1 Distribución T

2 Comparación de dos medias con t de Student

2.1 Supuestos de la prueba

2.2 Establecer las hipótesis

2.3 Determinar el tipo de test; una o dos colas

2.4 Formula

3 Ejemplo (manual):

4 Prueba T para una sola muestra

4.1 Ejercicio:

5 Comparación de dos medias independientes

5.1 Ejemplo 1: ¿Hay diferencia entre la longitud del culmen entre las poblaciones de pingünos?

5.2 Ejemplo 2: ¿ Hay diferencia en la longitud de la aleta entre pingünio Adelie y Gentoo?

5.3 Ejemplo 3: ¿Hay diferencia en el peso de los pokemon de roca y los de planta?

5.4 Ejemplo 4: ¿Hay diferencia entre la actividad enzimática de la citrato sintasa a 18°C y 32°C?; Prueba t con varianzas heterogeneas

6 Comparación de dos medias dependientes

6.1 Ejercicio: ¿Hay diferencia en el consumo de oxígeno (\(MO_2\)) tras la exposición a la temperatura?

7 Reportar los resultados de una prueba t con apa

7 Reportar los resultados de una prueba t con `apa`