Pruebas paramétricas y no paramétricas
Coriandrum sativum ha sido cultivada históricamente en nuestro país, pero los productores quieren saber si les convendría cambiar por el cultivo de C. tordylium. El INTA lleva a cabo un estudio a fin de comparar el rendimiento de ambas especies. De las áreas sembradas con cada una de ellas, se eligieron 15 parcelas al azar, se cosecharon en su totalidad y se determinó el rendimiento por cromatografía gaseosa (g/cm3).
En base a las observaciones realizadas, ¿podría el INTA recomendar a los productores cambiar de especie con un nivel de significación del 5%?
df_coriandrum <-read.delim2("Coriandrum.txt")
glimpse(df_coriandrum)Rows: 80
Columns: 2
$ Especie <chr> "C_sativum", "C_sativum", "C_sativum", "C_sativum", "C_sat~
$ Rendimiento <chr> "0.19", "0.201", "0.172", "0.193", "0.194", "0.225", "0.24~Cambiamos el tipo de dato de la variable Rendimiento
df_coriandrum$Rendimiento <- as.numeric(df_coriandrum$Rendimiento)Estadística Descriptiva
summary(df_coriandrum) Especie Rendimiento
Length:80 Min. :0.0410
Class :character 1st Qu.:0.1745
Mode :character Median :0.2120
Mean :0.2152
3rd Qu.:0.2525
Max. :0.3550 boxplot(Rendimiento~Especie,
data=df_coriandrum,
xlab="Especie",
ylab="Rendimiento",
col="orange",
border="brown"
)
Planteo de Hipótesis:
Ho:\(\mu_{ct} - \mu_{cs} = 0\)
Ha:\(\mu_{ct} - \mu_{cs}> 0\)
¿Conocemos las varianzas poblacionales?
Hipotesis para las varianzas poblacionales
Ho \(\frac{\sigma^2_{ct}}{\sigma^2_{cs}} =1\)
H1 \(\frac{\sigma^2_{ct}}{\sigma^2_{cs}} \neq 1\)
var.test(Rendimiento~Especie,
data=df_coriandrum,
alternative = "two.sided") # es una prueba bilateral
F test to compare two variances
data: Rendimiento by Especie
F = 0.86089, num df = 39, denom df = 39, p-value = 0.6423
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.4553221 1.6276938
sample estimates:
ratio of variances
0.8608862 Vuelvo a las Hipotesis de problema:
Ho:\(\mu_{ct} - \mu_{cs} = 0\)
Ha:\(\mu_{ct} - \mu_{cs}> 0\)
Con el parámetro var.equal se especifica si las varianzas de las poblaciones son similares o no, por ejemplo var.equal= T indica que las varianzas de ambas poblaciones son iguales.
En alternative se especifica el tipo de prueba de hipótesis.
t.test(Rendimiento~Especie ,
data=df_coriandrum,
alternative = "greater",
var.equal = T) #varianzas poblacionales iguales
Two Sample t-test
data: Rendimiento by Especie
t = -2.059, df = 78, p-value = 0.9786
alternative hypothesis: true difference in means between group C_sativum and group C_tordylium is greater than 0
95 percent confidence interval:
-0.05511309 Inf
sample estimates:
mean in group C_sativum mean in group C_tordylium
0.200000 0.230475 Si se desea obtener el límite inferior y superior del intervalo de confianza, debe correrse nuevamente la prueba, pero seleccionando la opción bilateral (alternative= two.sided)
¿Conclusiones?
En un trabajo de investigación se utilizaron 16 parcelas experimentales con dos plantas de avena cada una con el fin de estudiar el efecto promotor del crecimiento de una solución de potasio. En cada parcela, una planta elegida al azar fue tratada con la solución de potasio y la otra no (control).
df_avena <-read_excel("avena.xlsx")
glimpse(df_avena)Rows: 16
Columns: 3
$ parcelas <chr> "p1", "p2", "p3", "p4", "p5", "p6", "p7", "p8", "p9", "p10~
$ control <dbl> 22.37, 23.37, 28.65, 28.31, 23.53, 25.70, 26.46, 20.81, 25~
$ tratamiento <dbl> 21.46, 23.02, 21.18, 27.81, 26.78, 29.68, 25.32, 28.03, 26~\(\LARGE X_d=\) diferencia en el crecimiento de una avena tratada y una avena control
Hipotesis
Ho \(\mu_d = 0\)
Ha \(\mu_d > 0\)
df_avena <-df_avena %>%
mutate(diferencia= tratamiento-control)qqPlot(df_avena$diferencia)
[1] 3 16Prueba t para muestras pareadas
t.test(Pair(control, tratamiento) ~ 1,
alternative = "greater",
data = df_avena)
Paired t-test
data: Pair(control, tratamiento)
t = -0.50591, df = 15, p-value = 0.6899
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
-2.352575 Inf
sample estimates:
mean of the differences
-0.526875 df2_avena <-stack(df_avena[2:3]) |>
relocate (tratamiento=ind) |>
rename(altura=values)
head(df2_avena, 32) tratamiento altura
1 control 22.37
2 control 23.37
3 control 28.65
4 control 28.31
5 control 23.53
6 control 25.70
7 control 26.46
8 control 20.81
9 control 25.58
10 control 26.41
11 control 24.80
12 control 22.13
13 control 27.95
14 control 23.64
15 control 21.31
16 control 27.09
17 tratamiento 21.46
18 tratamiento 23.02
19 tratamiento 21.18
20 tratamiento 27.81
21 tratamiento 26.78
22 tratamiento 29.68
23 tratamiento 25.32
24 tratamiento 28.03
25 tratamiento 26.30
26 tratamiento 29.40
27 tratamiento 30.79
28 tratamiento 24.81
29 tratamiento 24.02
30 tratamiento 22.08
31 tratamiento 25.25
32 tratamiento 20.61t.test(df2_avena$altura~df2_avena$tratamiento,
alternative = "greater",
paired = T)
Paired t-test
data: df2_avena$altura by df2_avena$tratamiento
t = -0.50591, df = 15, p-value = 0.6899
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
-2.352575 Inf
sample estimates:
mean of the differences
-0.526875 En un estudio efectuado a fin de caracterizar la calidad y producción de aceite de oliva en la provincia de Catamarca de la República Argentina, se estudiaron dos de las variedades más conocidas. Para ello, se tomaron muestras de aceitunas de distintos ejemplares a una misma altura de copa de aproximadamente dos metros, y de todos los puntos cardinales de la misma, a efectos de evitar las variaciones debidas a la posición del fruto en la planta. Las aceitunas fueron secadas en estufa y se les determinó su contenido porcentual de aceite por extracción química.
# Cargamos los datos
Arbequina=c(34.5 , 20.1 , 21.8 , 18.2 , 19.5 , 20.2 , 22.5 , 23.9 , 22.1 , 24.2 )
Carolea=c(16.4 , 14.8 , 17.8 , 12.3 , 11.9 , 15.5 , 13.4 , 16, 15.8 , 16.2 )Realizamos la prueba de normalidad
shapiro.test(Arbequina) # Testea la normalidad de los datos
Shapiro-Wilk normality test
data: Arbequina
W = 0.76828, p-value = 0.00596shapiro.test(Carolea ) # Testea la normalidad de los datos
Shapiro-Wilk normality test
data: Carolea
W = 0.92481, p-value = 0.3989library(reshape2)
myList <-list(Arbequina, Carolea)
df <- melt(myList)
qplot(factor(L1), value, data = df, geom = "boxplot", xlab="Especies")
\(\theta_A\) a la mediana poblacional (posición central) del contenido de aceite de la variedad Arbequina.
\(\theta_C\) a la mediana poblacional (posición central) del contenido de aceite de la variedad Carolea.
Hipótesis
Ho \(\theta_A = \theta_C\)
Ha \(\theta_A \neq \theta_C\)
wilcox.test(Arbequina,
Carolea,
alternative = "two.sided")
Wilcoxon rank sum exact test
data: Arbequina and Carolea
W = 100, p-value = 1.083e-05
alternative hypothesis: true location shift is not equal to 0library(RVAideMemoire)# Paquete que contiene funciones misceláneas útiles en bioestadística
aceite<-read_excel("aceite.xlsx")
mood.medtest(Aceite~Variedad,
data=aceite)# Realiza el test de la mediana de Mood
Mood's median test
data: Aceite by Variedad
p-value = 1.083e-05Concluimos que rechazamos la hipótesis de la igualdad de las medianas de las dos variedades.
Interesa estudiar la capacidad detoxificadora del césped E. ophiuroides en suelos contaminados con Cd. A 20 macetas con césped se les asignó una de 5 dosis de Cd diferente (60, 120, 180, 240 y 300 mg Cd kg-1); 4 macetas por dosis. Se midió el Cd acumulado por la planta (expresado como mg Cd kg-1 materia seca)
#cargamos los datos
cadmio <- read.csv2("cadmio.csv")boxplot(cd_tallo~dosis_cd,
data=cadmio,
xlab="Dosis Cd (mg/kg suelo",
ylab="Cd tallo (mg/g MS)",
col="orange")
Hipótesis estadísticas
Ho \(\mu_{60} = \mu_{120} = \mu_{180}= \mu_{240}= \mu_{300}\)
Ha \(\mu_{i} \neq \mu\)
Las muestras deben ser aleatorias y las observaciones independientes
Las varianzas de las subpoblaciones deben ser iguales (homocedasticidad)
La distribución de cada subpoblación es normal
Definimos el modelo:
modelo<-aov(cd_tallo ~ factor(dosis_cd),
data=cadmio)
summary(modelo) Df Sum Sq Mean Sq F value Pr(>F)
factor(dosis_cd) 4 87303 21826 129 2e-11 ***
Residuals 15 2538 169
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Primero, se chequean los supuestos
Calculamos los residuos del modelo
e<-resid(modelo) # residuos
re<-rstandard(modelo) #residuos estandarizados
pre<-predict(modelo) #predichos
res<-cbind(cadmio$dosis_cd,cadmio$cd_tallo,pre,e,round(re,2))
colnames(res)<-c("dosis Cd", "Cd tallo", "Predichos", "Residuos", "residuos std")
res dosis Cd Cd tallo Predichos Residuos residuos std
1 60 23.2 30.30 -7.10 -0.63
2 60 16.2 30.30 -14.10 -1.25
3 60 52.7 30.30 22.40 1.99
4 60 29.1 30.30 -1.20 -0.11
5 120 52.5 54.75 -2.25 -0.20
6 120 45.7 54.75 -9.05 -0.80
7 120 52.9 54.75 -1.85 -0.16
8 120 67.9 54.75 13.15 1.17
9 180 133.5 140.00 -6.50 -0.58
10 180 156.9 140.00 16.90 1.50
11 180 123.9 140.00 -16.10 -1.43
12 180 145.7 140.00 5.70 0.51
13 240 166.8 168.50 -1.70 -0.15
14 240 165.9 168.50 -2.60 -0.23
15 240 184.3 168.50 15.80 1.40
16 240 157.0 168.50 -11.50 -1.02
17 300 208.4 202.30 6.10 0.54
18 300 189.9 202.30 -12.40 -1.10
19 300 217.7 202.30 15.40 1.37
20 300 193.2 202.30 -9.10 -0.81res<-as.data.frame(res)
par(mfrow = c(1, 2))
plot(pre, re, xlab="Predichos", ylab="Residuos estandarizados",main="Grafico de dispersion de RE vs PRED" )
abline(0,0)
qqPlot(e)
[1] 3 10Por prueba de hipótesis
shapiro.test(e)
Shapiro-Wilk normality test
data: e
W = 0.92721, p-value = 0.1364library(moments)
agostino.test(e)
D'Agostino skewness test
data: e
skew = 0.48189, z = 1.05937, p-value = 0.2894
alternative hypothesis: data have a skewnessPrueba de homogeneidad de varianzas
library(car)
cadmio$dosis_cd <- as.factor(cadmio$dosis_cd)
leveneTest(cd_tallo ~ dosis_cd,
data=cadmio)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 4 0.4306 0.7844
15 El test de Levene no rechaza la hipótesis nula de homocedasticidad.
#modelo<-aov(cd_tallo ~ factor(dosis_cd),
# data=cadmio)
summary(modelo) Df Sum Sq Mean Sq F value Pr(>F)
factor(dosis_cd) 4 87303 21826 129 2e-11 ***
Residuals 15 2538 169
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparaciones a posteriori
IC para la diferencia de dos medias: equivale a una PH para dif de medias (¿el cero pertenece al IC?).
TukeyHSD(modelo) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = cd_tallo ~ factor(dosis_cd), data = cadmio)
$`factor(dosis_cd)`
diff lwr upr p adj
120-60 24.45 -3.95090463 52.8509 0.1088118
180-60 109.70 81.29909537 138.1009 0.0000000
240-60 138.20 109.79909537 166.6009 0.0000000
300-60 172.00 143.59909537 200.4009 0.0000000
180-120 85.25 56.84909537 113.6509 0.0000012
240-120 113.75 85.34909537 142.1509 0.0000000
300-120 147.55 119.14909537 175.9509 0.0000000
240-180 28.50 0.09909537 56.9009 0.0490061
300-180 62.30 33.89909537 90.7009 0.0000533
300-240 33.80 5.39909537 62.2009 0.0163461plot(TukeyHSD(modelo))
Se realizó una intervención educativa innovadora para mejorar el rendimiento de los estudiantes. Dentro de los grupos de clasificación, el A es el grupo de control y los restantes, B y C, son los grupos con distintas innovaciones.
Puntajes=c(13,27,26,22,28,27,43,35,47,32,31,37,33,33,33,26,44,33,54)
Grupo=as.factor(c(rep("A",6), rep("B", 6), rep("C", 7)))
Rendimiento=data.frame (Grupo, Puntajes)library(pgirmess)
ggplot(Rendimiento,
aes(x=Grupo,y=Puntajes,fill=Grupo)) +
geom_boxplot() +
xlab("") +
scale_fill_brewer(palette="Pastel1")
# Produce boxplots
grupoA=Rendimiento[Rendimiento$Grupo=="A",2]
grupoB=Rendimiento[Rendimiento$Grupo=="B",2]
grupoC=Rendimiento[Rendimiento$Grupo=="C",2]shapiro.test(grupoA)
Shapiro-Wilk normality test
data: grupoA
W = 0.76163, p-value = 0.02583shapiro.test(grupoB)
Shapiro-Wilk normality test
data: grupoB
W = 0.92057, p-value = 0.5095shapiro.test(grupoC)
Shapiro-Wilk normality test
data: grupoC
W = 0.82769, p-value = 0.07607Ho Los tres grupos tienen la misma posición para la variable de estudio dada por el puntaje
Ha al menos un grupo tiene diferente posición para la variable en estudio dada por el puntaje.
kruskal.test(Puntajes, Grupo)# Realiza el test de Kruskal Wallis
Kruskal-Wallis rank sum test
data: Puntajes and Grupo
Kruskal-Wallis chi-squared = 9.9265, df = 2, p-value = 0.00699Comparación múltiple entre tratamientos luego del test de Kruskal Wallis
kruskalmc(Puntajes~Grupo )Multiple comparison test after Kruskal-Wallis
alpha: 0.05
Comparisons
obs.dif critical.dif stat.signif
A-B 9.250000 7.777876 TRUE
A-C 8.130952 7.494949 TRUE
B-C 1.119048 7.494949 FALSE