Pruebas de Hipótesis entre más de dos grupos - Análisis de Varianza (ANOVA)

En el tema anterior estuvimos haciendo comparaciones de medias (y varianzas) entre dos muestras; frecuentemente, sin embargo, necesitamos comparar las medias de más de dos muestras. La metodología estadística utilizada comúnmente para estos casos, es el análisis de varianza.

Modelo de Una-Vía (One-Way)

Tenemos el siguiente modelo para describir el valor de cualquier observación, proveniente de varios grupos (\(n_i, i=1...k\)), cada uno con \(n_{ij}\) observaciones:

\[y_{ij} = \mu + \alpha_i + e_{ij}\] donde:

\(\mu\) representa la media de todas las mediciones de todos los grupos.

\(\alpha_i\) representa la diferencia entre la media del grupo \(n_i\) y la media de todos los grupos.

\(e_{ij}\) representa el error aleatorio alrededor de los valores \(\mu + \alpha_i\), de una observación del grupo \(n_i\).

Hipótesis para el análisis de varianza de una vía

La hipótesis nula para el análisis de varianza de una vía, sería que las medias de los grupos son las mismas (e iguales a la media global, \(\mu\)), para esto todos los valores de \(\alpha_i\) deben ser igual a 0. La hipótesis alterna sería entonces, que al menos un \(\alpha_i\) debe ser diferente de 0.

Variaciones dentro de grupo y entre grupos, para la prueba de hipótesis

La siguiente ecuación es la base para la prueba de hipótesis que queremos:

\[y_{ij} - \bar{\bar y} = (y_{ij} - \bar y_i) + (\bar y_i - \bar{\bar y})\] donde:

\(y_{ij} - \bar{\bar y}\), representa la desviación de una medición con respecto a la media global.
\(y_{ij} - \bar y_i\), representa la desviación de una medición con respecto a la media de su grupo, y que se conoce como la variabilidad dentro de grupo.
\(\bar y_i - \bar{\bar y}\), representa las desviación de la media de un grupo con respecto a la media global, y es un indicador de la variabilidad entre grupos.

Las siguientes gráficas nos ilustran como los estimadores de la variabilidad nos permiten realizar la prueba de hipótesis:

(tomado de: Rosner, 2016)

En la figura anterior (a) la variabilidad entre grupos (B) es mucho mayor que dentro de grupos (A), por lo tanto uno o más \(\alpha_i\) son diferentes de 0, y la \(H_0\) podría ser rechazada.

(tomado de: Rosner, 2016)

En este caso (b), la variabilidad entre grupos (B) es muy baja, en comparación con A, y por lo tanto los valores de \(\alpha_i\) serían cercanos o iguales a 0, y por lo tanto no se puede rechazar la hipótesis nula.

Cálculo de la variabilidad en forma de Suma de Cuadrados (SS)

Para propósitos de cálculo del estadístico a ser usado para la prueba de hipótesis (F), se calculan los siguientes estimadores de la variabilidad total, dentro de grupos, y entre grupos:

Suma Total de Cuadrados (Total SS)

\[\sum_{i=1}^{k} \sum_{j=1}^{n_i}(y_{ij} - \bar{\bar y})^2\]

Suma de Cuadrados Dentro de Grupos (Within SS)

\[\sum_{i=1}^{k} \sum_{j=1}^{n_i}(y_{ij} - \bar y_i)^2\]

Suma de Cuadrados Entre Grupos (Between SS)

\[\sum_{i=1}^{k} \sum_{j=1}^{n_i}(\bar y_i - \bar{\bar y})^2\]

Cálculo del estadístico F para la prueba de hipótesis

La prueba de significancia para las hipótesis anteriores, se define como la razón entre la media de ‘Between SS’ (Between MS) y la media de ‘Within SS’ (Within MS).
Si esta razón es muy grande, se rechaza \(H_0\), y si es pequeña, no podemos rechazar \(H_0\).
La razón Between MS / Within MS se ajusta a la distribución F, con \(k-1\) y \(n-k\) grados de libertad.

Cálculo “rápido” de ‘Between SS’ y ‘Within SS’

\[Between\ SS = \sum_{i=1}^k n_i \bar y_{i}^2 - \frac{y..^2}{n}\] \[Within\ SS = \sum_{i=1}^k (n_i - 1)s_i^2\] donde:

\(y..\) es la suma de todas las observaciones de todos los grupos. \(k\), es el número de grupos. \(n\), el número total de observaciones de todos los grupos.
\(s_i^2\), la varianza ( de la muestra) de cada grupo.

Cálculo del estadístico F

Primero debemos calcular ‘Between MS’ y ‘Within MS’:

\[Between\ MS = \frac{Between\ SS}{k - 1}\]

\[Within\ MS = \frac{Within\ SS}{n - k}\]

y el estadístico F calculado:

\[F = \frac{Between\ MS}{Within\ MS}\]

Obtenemos el valor teórica: \(F_{k-1,n-k,1-\alpha}\)

Si \(F > F_{k-1,n-k,1-\alpha}\), rechazamos \(H_0\).
Si \(F < F_{k-1,n-k,1-\alpha}\), aceptamos \(H_0\).

(tomado de Rosner, 2016)

Ejemplo de tabla de resultados de prueba de hipótesis ANOVA de una vía

(tomado de Rosner, 2016)

Ejemplo 1

(tomado de Zar, 2014)
En una investigación, 19 cerdos jóvenes fueron asignados, al azar, a cuatro grupos experimentales. Cada grupo se alimentó con una dieta diferente (D1, D2, D3, D4). Luego de ser criados hasta adultos, se midió la masa corporal (kg) de cada animal. Queremos saber si la masa corporal resultó ser igual (\(H_0\)) para las cuatro dietas.

\[H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4\] \[H_A: la\ masa\ promedio\ de\ los\ cerdos\ no\ resultó\ igual\ en\ todas\ las\ dietas.\] \(\alpha = 0.05\)

Cálculos manuales con R

library(kableExtra)
#tamaño de muestra n_i
n1 <- length(D1)
n2 <- length(D2)
n3 <- length(D3)
n4 <- length(D4)
df.n <- data.frame(n1,n2,n3,n4)
kable(df.n, format = "markdown")

n1	n2	n3	n4
5	5	4	5

#sumatoria de grupos i
sumX1j <- sum(D1)
sumX2j <- sum(D2)
sumX3j <- sum(D3)
sumX4j <- sum(D4)
df.sumX <- data.frame(sumX1j,sumX2j,sumX3j,sumX4j)
kable(df.sumX, format = "markdown")

sumX1j	sumX2j	sumX3j	sumX4j
323.1	356.5	293.4	316.2

#media por grupo i
Xbar1 <- mean(D1)
Xbar2 <- mean(D2)
Xbar3 <- mean(D3)
Xbar4 <- mean(D4)
df.Xbar <- data.frame(Xbar1,Xbar2,Xbar3,Xbar4)
kable(df.Xbar, format = "markdown")

Xbar1	Xbar2	Xbar3	Xbar4
64.62	71.3	73.35	63.24

#suma de todas las mediciones
sumXij <- sum(df.sumX[1,])
#número total de mediciones
N <- sum(df.n[1,])
#media global
TotalXbar <- sumXij/N
df.Total <- data.frame(sumXij,N,TotalXbar)
kable(df.Total, format = "markdown")

sumXij	N	TotalXbar
1289.2	19	67.85263

#varianza de grupos i
varX1j <- var(D1)
varX2j <- var(D2)
varX3j <- var(D3)
varX4j <- var(D4)
df.varX <- data.frame(varX1j,varX2j,varX3j,varX4j)
#Total SS
TotalSS <- sum((cerditos - TotalXbar)^2, na.rm = TRUE)
#Group SS
GroupSS <- sum(df.n[1,]*(df.Xbar[1,]^2)) - (sumXij^2)/N
#Error SS
ErrorSS <- sum((df.n[1,] - 1)*df.varX[1,])
df.SS <- data.frame(TotalSS,GroupSS,ErrorSS)
kable(df.SS, format = "markdown")

TotalSS	GroupSS	ErrorSS
479.6874	338.9374	140.75

#grados de libertad y MS y F
v.groups <- length(df.n) - 1
v.error <- N - length(df.n)
GroupMS <- GroupSS/v.groups
ErrorMS <- ErrorSS/v.error
Fcalc <- GroupMS/ErrorMS
Falfa <- qf(0.95,v.groups,v.error)
Pr <- 1 - pf(Fcalc,v.groups,v.error)
df.prueba <- data.frame(v.groups,v.error,GroupMS,ErrorMS,Fcalc,Falfa,Pr)
kable(df.prueba, format = "markdown")

v.groups	v.error	GroupMS	ErrorMS	Fcalc	Falfa	Pr
3	15	112.9791	9.383333	12.0404	3.287382	0.000283

Usando función aov en R

#data frame con factor Dieta
df.cerdos <- data.frame(Dieta=c(rep("D1", times=length(D1)),
                        rep("D2", times=length(D2)),
                        rep("D3", times=length(D3)),
                        rep("D4", times=length(D4))),
                   masa=c(D1, D2, D3, D4))
analisis <- aov(masa ~ Dieta, df.cerdos)
summary(analisis)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Dieta        3  338.9  112.98   12.04 0.000283 ***
## Residuals   15  140.8    9.38                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Supuestos del ANOVA

La prueba de ANOVA, al igual que la prueba t, asume una distribución normal de los datos, y homocedasticidad (homogeneidad de varianzas). Por lo tanto, antes de continuar con otras pruebas es importante primero determinar si los datos cumplen con estos supuestos.

Prueba de normalidad

A continuación examinaremos gráficamente si los datos tienen una distribución normal, usando una gráfica de cuantiles-cuantiles (‘Q-Q plot’):

library(EnvStats)
qqPlot(df.cerdos$masa, add.line = TRUE, points.col = "blue", line.col = "red", )

A continuación la prueba de Shapiro-Wilk para normalidad, con una hipótesis nula de que los datos se ajustan a la distribución normal:

shapiro.test(df.cerdos$masa)

## 
##  Shapiro-Wilk normality test
## 
## data:  df.cerdos$masa
## W = 0.95808, p-value = 0.5352

No podemos rechazar la hipótesis nula con esta probabilidad (error tipo I), y aceptamos que los datos tienen una distribución normal.

Para la homogeneidad de varianza realizamos la prueba de Bartlett, que tiene una hipótesis nula de varianzas iguales:

bartlett.test(masa ~ Dieta, data=df.cerdos)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  masa by Dieta
## Bartlett's K-squared = 0.47515, df = 3, p-value = 0.9243

No podemos rechazar la hipótesis nula con esta probabilidad (error tipo I), y aceptamos que los datos por grupo tienen homogeneidad de varianza.

Pruebas post-hoc

El análisis de varianza nos indica, para cierto nivel de probabilidad, si un factor (o más de uno) produce desviaciones de un valor medio global; sin embargo no nos indica cuál o cuáles niveles del factor son los determinantes en la desviación.

Gráfica con intervalos de confianza

Mediante la siguiente gráfica, de medias e intervalos de confianza (95%), podemos hipotetizar cuáles dietas son diferentes entre sí y cuáles no.

library(plyr)
library(ggplot2)
ic.dietas <- ddply(df.cerdos, "Dieta", plyr::summarize,
                   masa.mean=mean(masa), masa.sd=sd(masa),
                   Lenght=NROW(masa),
                   tfrac=qt(p=.95, df=Lenght-1),
                   Lower=masa.mean - tfrac*masa.sd/sqrt(Lenght),
                   Upper=masa.mean + tfrac*masa.sd/sqrt(Lenght)
                   )
ggplot(ic.dietas, aes(x=masa.mean, y=Dieta)) + geom_point() + geom_errorbarh(aes(xmin=Lower, xmax=Upper), height=.3)

Para poder determinar si las diferencias específicas entre los niveles de los factores son significativas, debemos utilizar las llamadas pruebas post-hoc o procedimientos de comparaciones múltiples (MCP).

Utilizaremos la prueba de Tukey para probar cuáles son las dietas del Ejemplo 1 que se diferencian entre si.

library(agricolae)
posthoc <- TukeyHSD(analisis)
posthoc

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = masa ~ Dieta, data = df.cerdos)
## 
## $Dieta
##         diff        lwr       upr     p adj
## D2-D1   6.68   1.096263 12.263737 0.0168421
## D3-D1   8.73   2.807553 14.652447 0.0034914
## D4-D1  -1.38  -6.963737  4.203737 0.8906642
## D3-D2   2.05  -3.872447  7.972447 0.7530266
## D4-D2  -8.06 -13.643737 -2.476263 0.0041505
## D4-D3 -10.11 -16.032447 -4.187553 0.0009497

par(cex = 0.6)
plot(posthoc)

***

ANOVA con dos factores (‘two-way ANOVA’)

En el modelo anterior de ANOVA consideramos el efecto de un solo factor (dieta), vamos a estudiar el caso cuando queremos probar si dos factores son determinantes en las diferencias entre las medias de dos o más grupos (también llamados celdas).

Las hipótesis en este caso son más de una:

hipótesis acerca del efecto del factor A
hipótesis acerca del efecto del factor B
hipótesis acerca del efecto de la interacción entre los factores

Ejemplo 2

(tomado de Kabacoff, 2015)

Estaremos utilizando datos contenidos en el sistema base de R: ToothGrowth

TG <- ToothGrowth
head(TG)

tail(TG)

Conejillos de India (n = 60) fueron asignados al azar (10 por grupo, diseño balanceado) a tratamientos combinados (grupos) de tres dosis de ácido ascórbico (dose: 0.5, 1, o 2 mg) y dos tipos de bebidas (supp: jugo de naranja, OJ o bebida con vitamina C, VC). La variable respuesta, luego de un tiempo con los tratamientos, fue la longitud de los dientes frontales (len: mm). El modelo a usar es: \[len = interacción(supp\ x\ dosis)\]

library(gplots)  #para la gráfica

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

#tabla tratamientos
table(TG$supp,TG$dose)

##     
##      0.5  1  2
##   OJ  10 10 10
##   VC  10 10 10

#medias de los grupos
medias.grupos <- aggregate(TG$len, by=list(TG$supp,TG$dose), FUN=mean)
medias.grupos <- setNames(medias.grupos, c("Bebida", "Dosis, mg", "Media Long. dientes, mm"))
medias.grupos

#desviación estándar
sd.grupos <- aggregate(TG$len, by=list(TG$supp,TG$dose), FUN=sd)
sd.grupos <- setNames(sd.grupos, c("Bebida", "Dosis, mg", "DE Long. dientes, mm"))
sd.grupos

#convertir dosis (numérico) a factor
TG$dose <- factor(TG$dose)
#anova dos factores
anova2 <- aov(TG$len ~ TG$supp*TG$dose)
summary(anova2)

##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## TG$supp          1  205.4   205.4  15.572 0.000231 ***
## TG$dose          2 2426.4  1213.2  92.000  < 2e-16 ***
## TG$supp:TG$dose  2  108.3    54.2   4.107 0.021860 *  
## Residuals       54  712.1    13.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#gráfica de medias y intervalo confianza 95%
plotmeans(TG$len ~ interaction(TG$supp,TG$dose,sep = " "),
          connect = list(c(1,3,5),c(2,4,6)),
          col = c("orange", "darkgreen"),
          main = "Gráfica de Medias de Grupos con IC 95%",
          xlab = "Combinación de dosis y bebida",
          ylab = "Longitud de dientes frontales, mm")

La tabla de resultados del ANOVA indica que tanto los efectos principales (bebida y dosis) como la interacción entre ambos son significativos para las diferencias entre los grupos. La gráfica nos muestra los detalles de las diferencias entre grupos: a dosis altas no hay diferencias entre las bebidas, mientras que a dosis baja e intermedia si hay un mayor efecto con la bebida jugo de naranja.

ANOVA con Repeticiones

En este tipo de diseño, los mismos individuos son medidos más de una vez en el tiempo, o con un nuevo tratamiento.

Ejemplo 3

(tomado de Kabacoff, 2015)

Vamos a usar los datos de un experimento ecofisiológico, en el cual plantas (Plant) de la gramínea Echinochloa crus-galli, provenientes de dos lugares (Type: Quebec, Mississippi), fueron pretratadas (Treatment) con frío (‘chilled’) o no (‘non-chilled’). Luego las plantas fueron tratadas con diferentes concentraciones ambientales de dióxido de carbono (conc: 7 niveles, \(\mu L\ /\ L\)), y se midió su tasa de asimilación de C (uptake: \(\mu mol\ /\ m^2\ /\ seg.\)). Para el presente análisis solo se van a usar las plantas pretratadas (‘chilled’).

#selección de plantas enfriadas previamente`
fotosintesis <- subset(CO2, Treatment=='chilled')
fotosintesis

Para continuar el análisis debemos indicar que la concentración de dióxido de carbono (conc) es un factor. El modelo para la prueba es el siguiente: \[uptake = interacción(conc\ x\ Type) + error(Plant\ en\ conc)\]

#variable numérica (conc) a factor
fotosintesis$conc <- factor(fotosintesis$conc)
#anova con repeticiones
anovarep <- aov(uptake ~ conc*Type + Error(Plant/(conc)), fotosintesis)
summary(anovarep)

## 
## Error: Plant
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## Type       1 2667.2  2667.2   60.41 0.00148 **
## Residuals  4  176.6    44.1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Plant:conc
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## conc       6 1472.4  245.40   52.52 1.26e-12 ***
## conc:Type  6  428.8   71.47   15.30 3.75e-07 ***
## Residuals 24  112.1    4.67                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#gráfica de grupos
boxplot(uptake ~ Type*conc, data = fotosintesis, col = (c("blue", "red")),
        las=2, cex.axis=0.70,
        main = "Fotosíntesis de plantas de dos localidades",
        ylab = "Tasa de asimilación CO2, micromol/m^2/seg.",
        cex.lab=0.75)

Las tablas de resultados de ANOVA indican un efecto significativo del lugar de origen de las plantas (Type), e igualmente de la concentración, y su interacción con Type, en la tasa de asimilación de CO2 para estas plantas. La gráfica muestra que el efecto de la concentración de CO2 es mayor en las plantas que provienen de Quebec.

ANCOVA

Aunque en una prueba ANOVA de una vía estamos asumiendo que el único o mayor efecto sobre la variable respuesta, es el factor considerado o tratamiento, otras variables pueden actuar como pseudofactores, produciendo un efecto conjunto con el factor considerado, y las llamamos covariables. Vamos a examinar este tipo de análisis para un ANOVA de una vía (un factor), con una covariable. A este tipo de análisis se le conoce como análisis de covarianza (ANCOVA).

Ejemplo 4

(tomado de Kabacoff, 2015)

Ratonas preñadas fueron divididas en cuatro grupos de tratamientos, con diferentes dosis de una droga (dose: 0, 5, 50, o 500). El peso promedio de la camada (weigth) era la variable dependiente; el tiempo de gestación (gesttime) fue la covariable.

#paquete requerido con los datos
library(multcomp)
#datos litter y tabla de n por tratamientos
litter

table(litter$dose)

## 
##   0   5  50 500 
##  20  19  18  17

#media por grupo de tratamiento
medias.trat <- aggregate(litter$weight, by=list(litter$dose), FUN=mean)
medias.trat <- setNames(medias.trat, c("Dosis", "Masa, g"))
medias.trat

#ANCOVA
ancova <- aov(weight ~ gesttime + dose, data=litter)
summary(ancova)

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## gesttime     1  134.3  134.30   8.049 0.00597 **
## dose         3  137.1   45.71   2.739 0.04988 * 
## Residuals   69 1151.3   16.69                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los resultados muestran un efecto significativo del tiempo de gestación, y marginalmente significativo de la dosis. Esto quiere decir que las medias de masa calculadas para los grupos tratados, están afectadas por el tiempo de gestación; es posible remover este efecto, usando la función effect:

library(effects)
effect("dose", ancova)

## 
##  dose effect
## dose
##        0        5       50      500 
## 32.35367 28.87672 30.56614 29.33460

Pruebas no-paramétricas

Cuando el anális de varianza basado en el estadístico F no se puede aplicar (mediciones no tienen distribución normal y no hay homogeneidad de varianzas), la mejor alternativa no-paramétrica (basada en rangos) es la prueba de Kruskal-Wallis. En general, la prueba no determina si las medias (o medianas) son iguales, sino más bien si las distribuciones de donde las muestras fueron tomadas, son diferentes o no.

Los datos deben ordenarse dentro de cada grupo k (usualmente de menor a mayor), e indicar su rango de manera global (comparando todos los grupos). Si hay valores iguales, se asigna como rango el promedio de los rangos ocupados por los valores iguales (por ejemplo: si los primeros tres valores son iguales, se les asignará el valor de rango ligado (1+2+3)/3=2, a los tres valores iguales)

Para esta prueba se utiliza el estadístico \(H\), cuya fórmula es: \[H = \frac{12}{N(N+1)} \sum_{i=1}^k \frac{R_{i}^2}{n_i} - 3(N+1)\] Donde:

\(n_i\): número de obsevaciones en grupo i
\(N\): número total de observaciones en los k grupos
\(R_i\): suma de rangos en cada grupo i

Cuando hay valores con rangos ligados es necesario aplicar una corrección (\(C\)) al valor calculado de \(H\): \[C = 1 - \frac{\sum_{i=1}^m(t_{i}^3 - t_i)}{N^3 - N}\] Donde:

\(m\): es el número de grupos de rangos ligados
\(t_i\): es el número de valores en un rango ligado
\(N\): es número total de valores entre todos los grupos

La corrección sería: \[H_c = \frac{H}{C}\]

Los valores críticos para el estadístico \(H\) (o \(H_c\)) se obtienen de tablas para 3 o 4 grupos; para más grupos y muestras de tamaño grande se pueden utilizar los estadístico \(\chi^2\) o \(F\).

Para \(\chi^2\) se usa \(\nu = k-1\), y para F se utiliza la siguiente fórmula: \[F = \frac{(N-k)H_c}{(k-1)(N-1-H_c)}\] con \(\nu_1=k-1\ y\ \nu_2=N-k-1\)

Ejemplo 5

(tomado de Zar, 2014) Un limnólogo toma muestras en cuatro estanques, y mide el pH de cada muestra. Las hipótesis a probar son: \[H_0: el\ pH\ es\ igual\ en\ los\ cuatro\ estanques\] \[H_A: el\ pH\ no\ es\ el\ mismo\ en\ los\ cuatro\ estanques\] \(\alpha=0.05\)

#datos en data.frame
est1 <- c(7.68,7.69,7.7,7.7,7.72,7.73,7.73,7.76)
est2 <- c(7.71,7.73,7.74,7.74,7.78,7.78,7.8,7.81)
est3 <- c(7.74,7.75,7.77,7.78,7.8,7.81,7.84,"NA")
est4 <- c(7.71,7.71,7.74,7.79,7.81,7.85,7.87,7.91)
phdata <- data.frame(est1,est2,est3,est4)
phdata

Cálculos manuales

library(kableExtra)
#arreglo de los datos por factor estanque
#estanque 3 sin NA
est3na <- c(7.74,7.75,7.77,7.78,7.8,7.81,7.84)
phdata.est <- data.frame(Est=c(rep("est1", times=length(est1)),
                        rep("est2", times=length(est2)),
                        rep("est3", times=length(est3na)),
                        rep("est4", times=length(est4))),
                   pH=c(est1, est2, est3na, est4))
#asignación de rangos de menor a mayor, con promedio de iguales
RpH <- rank(phdata.est$pH, ties.method = "average")
#tabla con resultados
ranked.pH <- cbind(phdata.est,RpH)
ranked.pH

#suma rangos por estanque
Ri <- aggregate(ranked.pH$RpH, list(phdata.est$Est), sum)
kable(Ri, format = "markdown", col.names = c("Estanque", "Suma rangos"))

Estanque	Suma rangos
est1	55.0
est2	132.5
est3	145.0
est4	163.5

#tamaños de muestras y total
n1 <- length(est1)
n2 <- length(est2)
n3 <- length(est3na)
n4 <- length(est4)
N <- n1+n2+n3+n4
#cálculo de H
H <- ((12)/(N*(N+1)))*
  sum((Ri[1,2]^2)/n1 + (Ri[2,2]^2)/n2 + (Ri[3,2]^2)/n3 + (Ri[4,2]^2)/n4) - 3*(N+1)
#rangos ligados
rangos.ligados <- as.data.frame(table(ranked.pH$RpH))
rangos.ligados

#cálculo de sumatoria t
sum.t <- sum((2^3-2)+(3^3-3)+(3^3-3)+(4^3-4)+(3^3-3)+(2^3-2)+(3^3-3))
#cálculo de C
C <- 1 - sum.t/(N^3-N)
#calculo de Hc
Hc <- H/C
#chi-cuadrado
Chi2 <- qchisq(1-0.05,3)
#Fcalculado k:número de grupos
k = 4
Fcalc <- ((N-k)*Hc)/((k-1)*(N-1-Hc))
#Ftabla
Ftabla <- qf(1-0.05,k-1,N-k-1)
resulta <- data.frame(H,C,Hc,Chi2,Fcalc,Ftabla)
prueba <- kable(resulta, format = "markdown")
prueba

H	C	Hc	Chi2	Fcalc	Ftabla
11.87612	0.9943548	11.94354	7.814728	5.953097	2.975154

Podemos concluir que hay estanques con pH diferentes.

Usando función kruskal.test

phdata.est

kruskal.test(pH ~ Est, data = phdata.est)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  pH by Est
## Kruskal-Wallis chi-squared = 11.944, df = 3, p-value = 0.007579

Referencias

Kabacoff, R. 2015. R in Action. Data analysis and graphics with R, Second Ed. ed. Manning, Shelter Island.

Rosner, B. 2016. Fundamentals of biostatistics, 8th edition. ed. Cengage Learning, Boston, MA.

Zar, J.H. 2014. Biostatistical analysis. Pearson India, Noida.

ANOVA

D. S. Fernandez

10/14/2018

Pruebas de Hipótesis entre más de dos grupos - Análisis de Varianza (ANOVA)

Modelo de Una-Vía (One-Way)

Hipótesis para el análisis de varianza de una vía

Variaciones dentro de grupo y entre grupos, para la prueba de hipótesis

Cálculo de la variabilidad en forma de Suma de Cuadrados (SS)

Suma Total de Cuadrados (Total SS)

Suma de Cuadrados Dentro de Grupos (Within SS)

Suma de Cuadrados Entre Grupos (Between SS)

Cálculo del estadístico F para la prueba de hipótesis

Cálculo “rápido” de ‘Between SS’ y ‘Within SS’

Cálculo del estadístico F

Ejemplo 1

Cálculos manuales con R

Usando función aov en R

Supuestos del ANOVA

Prueba de normalidad

Pruebas post-hoc

Gráfica con intervalos de confianza

ANOVA con dos factores (‘two-way ANOVA’)

Ejemplo 2

ANOVA con Repeticiones

Ejemplo 3

ANCOVA

Ejemplo 4

Pruebas no-paramétricas

Ejemplo 5

Cálculos manuales

Usando función kruskal.test

Referencias