CLASE 8

INFERENCIA A PARTIR DE DOS MUESTRAS

1. OBJETIVO

En esta segunda fase del curso nos enfocaremos en formular conclusiones inferenciales de parámetros poblacionales a patir de datos muestrales. Usaremos los conceptos de intervalo de confianza y pruebas de hipótesis vistos en los capítulos anteriores.

2. PUNTOS CLAVE

Comparar entre proporciones de diferentes grupos

Parámetro	Muestra	Población
Proporción	\(\hat{p}\)	\(p\)
Media	\(\bar{x}\)	\(\mu\)
Desviación estandar	\(s\)	\(\sigma\)
Varianza	\(s^2\)	\(\sigma^2\)

	\(H_0\) es verdadera	\(H_0\) es falsa
Rechaza \(H_0\)	ERROR TIPO I (\(\alpha\))	Decisión correcta
Acepta \(H_0\)	Decisión correcta	ERROR TIPO II (\(\beta\))

3. TEMARIO

Dos proporciones.
Dos medias: muesras independientes
Dos muestras dependientes (pares relacionados)
Dos varianzas o desviaciones estándar

4. DOS PROPORCIONES

4.1. Notaciones para dos proporciones

Símbolo	Muestra	Población
Número de datos	\(n_1\) y \(n_2\)	\(N\)
Número de éxitos	\(x_1\) y \(x_2\)	\(m\)
Proporción	\(\hat{p}_1\) y \(\hat{p}_2\)	\(p_1\)
Complento de proporción	\(\hat{q}_1\) y \(\hat{q}_2\)	\(q_1\)
Proporcional muestral agrupada	\(\bar{p}=\frac{x_1 + x_2}{n_1 + n_2}\)	NA
Complemento proporcional muestral agrupada	\(\bar{q}= 1 - \bar{p}\)	NA

4.2 Requisitos o supuestos

Proporciones muestrales de muestras aleatorias
Muestras independientes
Para cada muestras hay por lo menos 5 éxitos y 5 fracasos

4.3 Estimación de valores críticos

Parámetro	Símbolo	Detalle
Margen de error	\(E\)	\(E=z_{\alpha/2}\sqrt{\frac{\hat{p}_1.\hat{q}_1}{n_1}+\frac{\hat{p}_2.\hat{q}_2}{n_2}}\)
Porporción	\(\hat{p}\)	\(\hat{p} = \frac{x_1}{n_1}\)
Tipo de distribución	\(z\)	Distribución normal estandar

4.4 Ejemplo

Ahora vamos a trabajar con la base de datos de la clase anterior de “alevin”, para lo cual lo tenemos que importar de nuestra carpeta

alevin <- read.delim("C:/Users/acard/CloudStation/R-project/bioestadistica_2020/chp8/data/alevin.txt")
head(alevin)

##   REPETICION POZO ESTADIO       SUBSTANCIA CONCENTRACION TIEMPO MORTALIDAD
## 1          1    A  ALEVIN casiopeina_IIIia             0      0          0
## 2          1    B  ALEVIN casiopeina_IIIia             0      0          0
## 3          1    A  ALEVIN casiopeina_IIIia             1      0          0
## 4          1    B  ALEVIN casiopeina_IIIia             1      0          0
## 5          1    A  ALEVIN casiopeina_IIIia            10      0          0
## 6          1    B  ALEVIN casiopeina_IIIia            10      0          0
##   MALFORMACION MOVILIDAD TOTAL
## 1            0         5     5
## 2            0         5     5
## 3            0         5     5
## 4            0         5     5
## 5            0         5     5
## 6            0         5     5

Ahora calculemos si la porporciónn de mortalidad de alevines es la misma a las 72 horas de ser expuestos a “casiopeina_IIIia” a concentraciones mayores o menores a 20 ng/mL. Lo primero que tenemos que hacer es generar nuestro subset de datos

h.72 <- subset(x = alevin, 
               subset = TIEMPO == "72" & SUBSTANCIA == "casiopeina_IIIia")
head(h.72)

##    REPETICION POZO ESTADIO       SUBSTANCIA CONCENTRACION TIEMPO MORTALIDAD
## 43          1    A  ALEVIN casiopeina_IIIia             0     72          1
## 44          1    B  ALEVIN casiopeina_IIIia             0     72          0
## 45          1    A  ALEVIN casiopeina_IIIia             1     72          1
## 46          1    B  ALEVIN casiopeina_IIIia             1     72          0
## 47          1    A  ALEVIN casiopeina_IIIia            10     72          0
## 48          1    B  ALEVIN casiopeina_IIIia            10     72          0
##    MALFORMACION MOVILIDAD TOTAL
## 43            0         4     5
## 44            0         4     5
## 45            0         4     5
## 46            0         5     5
## 47            0         5     5
## 48            3         5     5

table(h.72$MORTALIDAD, h.72$CONCENTRACION)

##    
##     0 1 10 25 50 100 200
##   0 3 3  3  0  0   0   0
##   1 3 2  0  0  0   0   0
##   2 0 1  2  0  0   0   0
##   3 0 0  1  0  0   0   0
##   4 0 0  0  2  0   0   0
##   5 0 0  0  4  6   6   6

Paso #1

Ahora debemos revisar son los requisitos de la prueba: + Los datos provienen de muestras aleatorias. + Los datos son independientes + El número de exitos y fracasos es mayor a 5

Paso #2

Establecemos la prueba de hipótesis \[H_0:\hat{p_1}=\hat{p_2}\] \[H_1:\hat{p_1}\neq \hat{p_2}\] #### Paso #3 Estimamos el valor de la porporción muestral acumulada que se define por la siguiente fórmula \[\hat{p} = \frac{x_1+x_2}{n_1+n_2}\] Para obtner la suma de los eventos exitosis \(x_i\) y el número total de animales \(n_i\) podemos usar la función tapply(). Por ejemplo, para obtener la suma de los eventos exitosos por concentración, usaríamos la siguiente función

x <- tapply(X = h.72$MORTALIDAD, INDEX = h.72$CONCENTRACION, FUN = sum)
x

##   0   1  10  25  50 100 200 
##   3   4   7  28  30  30  30

Ahora derminamos \(x_1\) y \(x_2\) si tenemos como punto de corte los 20ng/mL

x[1:3]

##  0  1 10 
##  3  4  7

x1 <- sum(x[1:3])
x2 <- sum(x[4:length(x)])

Ahora calculemos la n para cada una de las dos muestras

n <- tapply(X = h.72$TOTAL, INDEX = h.72$CONCENTRACION, FUN = sum)
n1 <- sum(n[1:3])
n2 <- sum(n[4:length(n)])

Ahora calculemos la porporción muestral acumulada

p.ac <- (x1+x2)/(n1+n2)
p.ac

## [1] 0.6285714

p.com <- 1-p.ac
p.com

## [1] 0.3714286

Paso #4

Ahora determinemos el valor de z con la fórmula \[z = \frac{(\hat{p}_1-\hat{p}_2)-(p_1-p_2)}{\sqrt{\frac{\bar{p}\bar{q}}{n_1}+\frac{\bar{p}\bar{q}}{n_2}}}\]

En la hipótesis nula asumimos que \(p_1=p_2\). Reemplazando los valores en la fórmula tenemos que

z <- ((x1/n1-x2/n2)-(0))/sqrt(((p.ac*p.com)/n1)+((p.ac*p.com)/n2))
z

## [1] -12.28573

El área bajo la curva para este valor de z lo podemos calcular con la función pnorm

pnorm(z)

## [1] 5.403152e-35

Paso #5

Determinmos los valores críticos del intervalo de confianza en la distribución z y determinamos si el valor del estadistico está en la zona de rechazo. Para calcular los valores críticos con un intervalo de confianza de 95% a dos colas, podemos utilizar la función qnorm()

# cola derecha
qnorm(0.025)

## [1] -1.959964

# cola izquierda
qnorm(0.025, lower.tail = F)

## [1] 1.959964

Entonces los valores críticos para el intevalo de confianza de 95% son (-1.96,1.96), es decir que el área de rechazo de la hipótesis nula está por debajo del valor z de -1.96 y por encima del valor de z de 1.96. El valor del estadistico de prueba es de -12.28, es decir, que está en la zona de rechazo de la hipótesis nula. Hagamos una gráfica de barras para evidenciar la diferencia

center <- barplot(c(x1/n1,x2/n2), main = "Proporción de mortalidad a 72h", col = colores1, las = 1, ylim = c(0,1), ylab = "Proporción", xlab = "Grupos")
axis(side = 1, at = center, labels = c("<20ng/mL", ">20ng/mL"))

R nos permite hacer este análisis con una función prop.test la cual tiene 4 parámetros que debemos asignar: + x= es el número de éxitos que tenemos en cada muestra + n= es el número de ensayos que se realizaron en cada muestra + conf.level= es el nivel de significancia de la prueba + alternative= indicamos si el área de rechazo es a dos colas, o una cola. Las opciones son: * “two.side” = a dos colas * “greater” = mayor que * “less” = menor que

prop.test(x = c(x1,x2), n = c(n1,n2),conf.level = 0.95,alternative = "two.sided")

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(x1, x2) out of c(n1, n2)
## X-squared = 147.41, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.9158031 -0.7397524
## sample estimates:
##    prop 1    prop 2 
## 0.1555556 0.9833333

R tambien nos permite generar una tabla de porporciones de forma sensilla con la función prop.table(). Para generar esta tabla, necesitamos ingresar una tabla con las frecuencias abosultas de los eventos de exito y fracaso ( en este caso, sería muertos y vivos respectivamente). Esta tabla de frecuencias absolutas la podemos generar con la función table() o con ayuda de la función tapply(). Para este ejemplo usaremos esta ultima. Para calcular el número de muertos por concentración (como ya lo habíamos hecho anteriormente) usamos la MORTALIDAD en el parámetro de “X”, y como factor de agrupación CONCENTRACIÓN en el parámetro “INDEX”, finalmente le indicamos a la función que sume los valores así agrupados indicando “sum” en el parámetro FUN. Para calcular el número de vivos, lo que hacemos es restar el número de muertos del total de individuos. Finalmente generamos una matriz de estos datos con la función cbind()

muertos <- tapply(X = h.72$MORTALIDAD, INDEX = h.72$CONCENTRACION, FUN = sum)
vivos <- tapply(X = h.72$TOTAL, INDEX = h.72$CONCENTRACION, FUN = sum) - muertos
db <- cbind(muertos, vivos)
db

##     muertos vivos
## 0         3    27
## 1         4    26
## 10        7    23
## 25       28     2
## 50       30     0
## 100      30     0
## 200      30     0

Ahora que ya tenemos nuestra tabla de frecuenias absolutas, podemos convertirlo a porporciones o frecuencias relativas usando la función prop.table(). Esta función necesita que definamos dos parámetros: + x= donde se denfine la tabla de frecuencias absolutas antes generada + margin= donde se define si los calculos se hacen horizontales (1) o verticales (2)

prop.table(x = db, margin = 1)

##       muertos      vivos
## 0   0.1000000 0.90000000
## 1   0.1333333 0.86666667
## 10  0.2333333 0.76666667
## 25  0.9333333 0.06666667
## 50  1.0000000 0.00000000
## 100 1.0000000 0.00000000
## 200 1.0000000 0.00000000

De la misma forma, podemos generar una tabla de porporciones con la función proportions(). En la cual definimos los parametros al igual que la anterior.

proportions(x = db,margin = 1)

##       muertos      vivos
## 0   0.1000000 0.90000000
## 1   0.1333333 0.86666667
## 10  0.2333333 0.76666667
## 25  0.9333333 0.06666667
## 50  1.0000000 0.00000000
## 100 1.0000000 0.00000000
## 200 1.0000000 0.00000000

5. DOS MEDIAS: MUESTRAS INDEPENDIENTES

5.1 Notación

Parámetro	Poblacional	Muestral
Media	\(\mu_1\)	\(\bar{x}_1\)
Desviación estandar	\(\sigma_1\)	\(s_1\)
Tamaño	\(N\)	\(n_1\)

5.2 Requisitos

Desviaciones estandar de las poblaciones desconocidas y no asumimos que son iguales.
Muestras independientes, aleatorias y simples
La n es > 30 o la población de donde procede la muestra tiene una distribución normal

5.3 Pueba estadistica para comparar dos medias

\[t=\frac{(\bar{x}_1 - \bar{x}_2)-(\mu_1-\mu_2)}{\sqrt{\frac{\sigma^2_1}{m_1}+\frac{\sigma^2_2}{n_2}}}\]

5.4 Ejemplo

En este caso, vamos a usar R, para determinar la diferencia entre medias. R tiene una función llamada t.test(), y la base de datos de R “iris”.

data("iris")
head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

Esta base de datos tiene 3 variables

levels(iris$Species)

## [1] "setosa"     "versicolor" "virginica"

Para hacer nuestra prueba vamos a hacer un subset de setosa y versicolor, usando la función subset()

prueba <- subset(x = iris, subset = Species == "setosa" | Species == "versicolor")
prueba$Species <- droplevels(prueba$Species)
boxplot(prueba$Sepal.Length~prueba$Species)

t.test(formula =prueba$Sepal.Length~prueba$Species)

## 
##  Welch Two Sample t-test
## 
## data:  prueba$Sepal.Length by prueba$Species
## t = -10.521, df = 86.538, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.1057074 -0.7542926
## sample estimates:
##     mean in group setosa mean in group versicolor 
##                    5.006                    5.936

Que es lo mismo si introducimos los valores de ambas especies por separado

t.test(x = prueba$Sepal.Length[prueba$Species == "setosa"],y = prueba$Sepal.Length[prueba$Species == "versicolor"])

## 
##  Welch Two Sample t-test
## 
## data:  prueba$Sepal.Length[prueba$Species == "setosa"] and prueba$Sepal.Length[prueba$Species == "versicolor"]
## t = -10.521, df = 86.538, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.1057074 -0.7542926
## sample estimates:
## mean of x mean of y 
##     5.006     5.936

La función t.test() tiene varios parámetros que podemos definir:

x= vector numérico de los valores de un grupo
y= vector numérico de los valores del otro grupo
formula= formula donde se define el vector númerico con el factor que las agrupa, se escribe de así: “data\(valor ~ data\)factor”
alternative= aquí se define el area de rechazo
- two.side dos colas
- greater cola derecha
- less cola izquierda
paired= aquí se define si los datos son pareados o no
conf.level= aquí se define el intervalo de confianza

Si reescribimos las funciones anteriores, para definir todos estos parámetros

t.test(x = iris$Sepal.Length[iris$Species=="setosa"],
       y = iris$Sepal.Length[iris$Species=="versicolor"], 
       alternative = "two.side", 
       paired = FALSE, 
       conf.level = 0.95)

## 
##  Welch Two Sample t-test
## 
## data:  iris$Sepal.Length[iris$Species == "setosa"] and iris$Sepal.Length[iris$Species == "versicolor"]
## t = -10.521, df = 86.538, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.1057074 -0.7542926
## sample estimates:
## mean of x mean of y 
##     5.006     5.936

6. DOS MUESTRAS DEPENDIENTES

En un buen diseño experimental, podemos definir facilmente los grupos pareados. Si por ejemplo, queremos evaluar la eficiencia de nuestras clases, podríamos tomar un exámen, antes y despues de la clase. Podríamos comparar la media de las notas antes y despues para ver si se asimiló los conocimientos adecuadamente. Es obvio ver que un mismo alumno tendrá dos notas, una antes y otra depues de la clase. Esto corresponde a grupos pareados.

6.1 Notación

Parámetro	Poblacional	Muestral
Diferencia individual	\(NA\)	\(d\)
Diferencia grupal	\(\mu_d\)	\(\bar{d}\)
Desviación estandar para las diferencias	\(\sigma_d\)	\(S_d\)
Número de pares de datos	\(N\)	\(n\)

6.2 Requisistos

Datos muestrales dependientes.
Muestras aleatorias simples
Número de diferencicas mayor a 30 o las diferencias tienen una distribución normal

6.3 Estadistico de prueba

\[t=\frac{\bar{d}-\mu_d}{\frac{S_d}{\sqrt{n}}}\]

6.4 Ejemplo

Para este ejemplo vamos a utilizar la base de datos de crecimiento bacteriano. Vamos a comparar el crecimiento de la especie “Chloridium” en los tiempos T2 y T5 , cuando tienen un ciclo de luz/oscuridad de 12h. Lo primero que tenemos que hacer es cargar nuesta base de datos

growth <- read.delim("C:/Users/acard/CloudStation/R-project/bioestadistica_2020/chp4/data/curva de crecimiento final.txt", stringsAsFactors=TRUE)
head(growth)

##   tiempo1   reactor tiempo2         sp Medio.BBM Medio.BG.11 Medio.F.2
## 1     T:0 reactor 1     12h Chloridium   819.222     554.778   554.778
## 2     T:1 reactor 1     12h Chloridium  1567.000    1248.222  1393.667
## 3     T:2 reactor 1     12h Chloridium  2093.667    3505.889  2392.556
## 4     T:3 reactor 1     12h Chloridium  2585.889    1912.556  1965.889
## 5     T:4 reactor 1     12h Chloridium  4985.111    5916.222 17935.000
## 6     T:5 reactor 1     12h Chloridium  5217.667    1925.889 11873.889

Ahora agamos un subset de nuestros datos de interes con la función subset()

df <- subset(x = growth,subset = growth$sp == "Chloridium" & growth$tiempo2 == "12h", select = c("tiempo1","reactor","Medio.BG.11") )
head(df,10)

##    tiempo1   reactor Medio.BG.11
## 1      T:0 reactor 1     554.778
## 2      T:1 reactor 1    1248.222
## 3      T:2 reactor 1    3505.889
## 4      T:3 reactor 1    1912.556
## 5      T:4 reactor 1    5916.222
## 6      T:5 reactor 1    1925.889
## 7      T:6 reactor 1    8968.333
## 8      T:0 reactor 2     554.778
## 9      T:1 reactor 2    3220.333
## 10     T:2 reactor 2    5934.000

Cambiamos el formato de la base de datos

df <- reshape(data = df, timevar = "tiempo1", idvar = "reactor", direction = "wide")
names(df) <- c("reactor","T0","T1","T2","T3","T4","T5","T6")
df

##      reactor       T0       T1       T2       T3        T4       T5        T6
## 1  reactor 1  554.778 1248.222 3505.889 1912.556  5916.222 1925.889  8968.333
## 8  reactor 2  554.778 3220.333 5934.000 5834.333 17168.667 7051.667  8602.000
## 15 reactor 3 2232.556 1661.444 1584.778 1681.444  2172.556 2953.667 11412.778

Ahora solo seleccionemos los datos que nos interesan

df<- df[,c("reactor","T2","T5")]
df

##      reactor       T2       T5
## 1  reactor 1 3505.889 1925.889
## 8  reactor 2 5934.000 7051.667
## 15 reactor 3 1584.778 2953.667

Para determinar la diferencia entre las muestras \(d\) lo cálculamos facilmente

df$d <- df$T5-df$T2
df

##      reactor       T2       T5         d
## 1  reactor 1 3505.889 1925.889 -1580.000
## 8  reactor 2 5934.000 7051.667  1117.667
## 15 reactor 3 1584.778 2953.667  1368.889

Para determinar la normalidad de la distribución de las diferencias de los datos podemos utilizar la prueba de normalidad de wilcoxon

wilcox.test(df$d)

## 
##  Wilcoxon signed rank exact test
## 
## data:  df$d
## V = 3, p-value = 1
## alternative hypothesis: true location is not equal to 0

chapter_9

Alfredo Cárdenas-Rivera

12/10/2020