Se seleccionaron muestras aleatorias de cada uno de tres tipos de individuos: empleados de una empresa fabricante de plaguicidas con un alto grado de exposición, agricultores expuestos al plaguicida por varias semanas cada año y personas sin exposición conocida al plaguicida. Se hicieron determinaciones de la acetilcolinesterasa en muestras de sangre de cada persona, obteniéndose los siguientes resultados.
¿Proporcionan estos datos evidencia suficiente que indique una diferencia, en promedio, entre los tres grupos de personas?
empleados<-c(6.4,6.6,6.8,6.9,9.5,6.1,7.5,8.2,4.1,5.5)
agricultores<-c(6.5,6.8,7,7.1,9.7,6.2,7.7,8.4,4.2,5.6)
noexp<-c(7.3,7.5,7.8,7.9,10.8,6.9,8.5,9.4,4.6,6.3)
length(empleados)[1] 10length(agricultores)[1] 10length(noexp)[1] 10grupo <- as.factor(rep(c("empleados", "agricultores", "noexp"), each = 10))
niveles <- c(empleados, agricultores, noexp)
length(grupo)[1] 30length(niveles)[1] 30data_leven<-data.frame(grupo,niveles)
data<-data.frame(empleados,agricultores,noexp)# Gráfico de caja y bigotes
boxplot(data, main = "Distribución de Acetilcolinesterasa en Sangre",names = c("Empleados", "Agricultores", "No Expuestos"), xlab = "Grupos de muestra",ylab = "Nivel en Sangre", ylim = c(4,11), col = c("darkslategray1","darkslategray3","darkslategray"))
# Histograma para cada grupo
par(mfrow = c(1,3))
hist(data$empleados, main = "Empleados", xlab = "Nivel de acetilcolinesterasa",ylab="Frecuencia", col = "darkslategray1", xlim =c(4,11))
hist(data$agricultores, main = "Agricultores", xlab = "Nivel de acetilcolinesterasa", ylab="Frecuencia", col = "darkslategray3", ylim =c(0,5),xlim =c(4,11))
hist(data$noexp, main = "No Expuestos", xlab = "Nivel de acetilcolinesterasa",ylab="Frecuencia", col = "darkslategray",ylim =c(0,5),xlim =c(4,11))
par(mfrow = c(1,1))summary(data) empleados agricultores noexp
Min. :4.100 Min. :4.200 Min. : 4.60
1st Qu.:6.175 1st Qu.:6.275 1st Qu.: 7.00
Median :6.700 Median :6.900 Median : 7.65
Mean :6.760 Mean :6.920 Mean : 7.70
3rd Qu.:7.350 3rd Qu.:7.550 3rd Qu.: 8.35
Max. :9.500 Max. :9.700 Max. :10.80 No hay diferencias en los niveles de acetilcolinesterasa promedio entre los tres grupos.
Al menos un grupo tiene una media diferente en los niveles de acetilcolinesterasa.
#la muestra tiene tamaño menr que 30 entonces usaremos shapiro test para comprobar normalidad.
shapiro.test(data$empleados) # p>0.05
Shapiro-Wilk normality test
data: data$empleados
W = 0.97778, p-value = 0.9522shapiro.test(data$agricultores)
Shapiro-Wilk normality test
data: data$agricultores
W = 0.98297, p-value = 0.9791shapiro.test(data$noexp)
Shapiro-Wilk normality test
data: data$noexp
W = 0.97871, p-value = 0.9579# con una significancia del 0.05 podemos afirmar que los datos son normales.library(car)Loading required package: carDataleveneTest(niveles ~ grupo, data=data_leven)#p>0.05Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 2 0.0626 0.9395
27 # Hay homogeneidad de varianzas con una significancia de 0.05Las muestras son independientes pues la medición se realizo en grupos distintos.
media_global <- mean(data_leven$niveles)
SST <- sum((niveles - media_global)^2)
SST<-round(SST,2)media_empleados <- mean(empleados)
media_agricultores <- mean(agricultores)
media_noexp <- mean(noexp)
medias <- c(media_empleados,media_agricultores,media_noexp)
# Suma de Cuadrados entre Grupos (SSA)
SSA<-sum(10*(medias-media_global)^2) # 10 porque son igual el tamaño de muestra para los tres grupos
SSA<-round(SSA,2)SSE <- sum((data$empleados - media_empleados)^2) +
sum((data$agricultores - media_agricultores)^2) +
sum((data$noexp - media_noexp)^2)
SSE<-round(SSE,2)SST[1] 70.48SSA[1] 5.06SSE[1] 65.42SST == SSA + SSE #se debe cumplir esto[1] TRUEk <- 3 # Número de grupos
N <- 30 # Tamaño total de la muestra 10 cada uno
# Cuadrado Medio entre Grupos (MSA)
MSA <- round(SSA / (k - 1),2)
MSA[1] 2.53# Cuadrado Medio dentro de los Grupos (MSE)
MSE <- round(SSE / (N - k),2)
MSE[1] 2.42# Estadístico F
rF <- round(MSA / MSE,2)
rF[1] 1.05Fuente_de_variación<-c("Tratamientos(entre grupos)","Error(dentro de grupos)","Total")
Suma_de_cuadrados<-c(SSA,SSE,SST)
Grados_de_libertad<-c(k-1,N-k,N-1)
Cuadrado_medio<-c(MSA,MSE,"")
RazónF<-c(rF,"","")
data.frame(Fuente_de_variación,Suma_de_cuadrados,Grados_de_libertad,Cuadrado_medio,RazónF) Fuente_de_variación Suma_de_cuadrados Grados_de_libertad
1 Tratamientos(entre grupos) 5.06 2
2 Error(dentro de grupos) 65.42 27
3 Total 70.48 29
Cuadrado_medio RazónF
1 2.53 1.05
2 2.42
3 # Funciones para cálculos ANOVA
calcular_SST <- function(data, niveles) {
media_global <- mean(data[[niveles]])
SST <- sum((data[[niveles]] - media_global)^2)
return(round(SST, 2))
}
calcular_SSA <- function(data, grupos, niveles) {
medias <- tapply(data[[niveles]], data[[grupos]], mean) # Medias por grupo
media_global <- mean(data[[niveles]]) # Media global
# Obtiene el tamaño de cada grupo
tam_grupos <- table(data[[grupos]])
# Calcula la SSA correctamente
SSA <- sum(tam_grupos * (medias - media_global)^2)
return(round(SSA, 2))
}
calcular_SSE <- function(data, grupos, niveles) {
medias <- tapply(data[[niveles]], data[[grupos]], mean)
SSE <- sum(sapply(unique(data[[grupos]]), function(g) {
sum((data[data[[grupos]] == g, niveles] - medias[g])^2)
}))
return(round(SSE, 2))
}
calcular_MSA <- function(SSA, k) {
MSA <- SSA / (k - 1)
return(round(MSA, 2))
}
calcular_MSE <- function(SSE, N, k) {
MSE <- SSE / (N - k)
return(round(MSE, 2))
}
calcular_F <- function(MSA, MSE) {
rF <- MSA / MSE
return(round(rF, 2))
}
generar_tabla_ANOVA <- function(data, grupos, niveles) {
k <- length(unique(data[[grupos]])) # Número de grupos
N <- nrow(data) # Tamaño total de la muestra
SST <- calcular_SST(data, niveles)
SSA <- calcular_SSA(data, grupos, niveles)
SSE <- calcular_SSE(data, grupos, niveles)
MSA <- calcular_MSA(SSA, k)
MSE <- calcular_MSE(SSE, N, k)
rF <- calcular_F(MSA, MSE)
Fuente_de_variación <- c("Tratamientos(entre grupos)", "Error(dentro de grupos)", "Total")
Suma_de_cuadrados <- c(SSA, SSE, SST)
Grados_de_libertad <- c(k-1, N-k, N-1)
Cuadrado_medio <- c(MSA, MSE, "")
RazónF <- c(rF, "", "")
tabla_anova <- data.frame(Fuente_de_variación, Suma_de_cuadrados, Grados_de_libertad, Cuadrado_medio, RazónF)
return(tabla_anova)
}
# Uso de la función
tabla_anova <- generar_tabla_ANOVA(data_leven, "grupo", "niveles")
print(tabla_anova) Fuente_de_variación Suma_de_cuadrados Grados_de_libertad
1 Tratamientos(entre grupos) 5.06 2
2 Error(dentro de grupos) 65.42 27
3 Total 70.48 29
Cuadrado_medio RazónF
1 2.53 1.05
2 2.42
3 anova<-aov(niveles~grupo,data=data_leven)
print(anova)Call:
aov(formula = niveles ~ grupo, data = data_leven)
Terms:
grupo Residuals
Sum of Squares 5.05867 65.42000
Deg. of Freedom 2 27
Residual standard error: 1.556587
Estimated effects may be unbalanced