PARCIAL DISEÑO DE EXPERIMENTOS
INTEGRANTES:
GINNA VALENTINA ROMERO RINCON 1026596100 (gvromeror@unal.edu.co)
LAURA CAMILA ALVAREZ GACHA 1026303112 (laucalvarezgac@unal.edu.co)
Punto 1
a.1)
En un estudio conducido en un ambiente controlado se tuvieron 72 macetas, cada una con una planta a a que a cierta edad se le midió el contenido de clorofila (indice de clorofila) con un sensor (SPAD). El total de macetas se correspondio con 9 tratamientos asociados a estrés hídrico. Se sabe que la varianza de las 72 observaciones es 823 y que la suma de cuadrados de los tratamientos (SCtrt) es 6000 Con esta información complete la tabla del ANOVA.
solución
\[¿Qué~ datos ~tenemos~ y ~ podemos ~ deducir?\]
\[n=72\\ tratamientos~~(trt)=9\\ repeticiones ~~(rpp)= \frac{n}{trt}= \frac{72}{9}= 8 \\ Varianza~~ (S^2)=823\\ SC_{trt}=6000\\ Grados~ de ~libertad~ de ~los~ tratamientos~~ (gl_{trt})=trt-1= 9-1=8\\ Grados~de~libertad~del~error~ (gl_{error})= trt*(rpp-1)=9(8-1)=63\\ Grados ~de ~libertad~ totales ~~ (gl_{totales})= n-1=72-1=71\]
n=72
trt=9
rpp=8
S2=823
SCtrt=6000
gl_trt=8
gl_error=63
gl_totales= 71para poder completar el ANOVA nos hace falta la suma de cuadrado de las repeticiones o suma de cuadrados del error (SCerror), y para obtener este valor partimos de esta fórmula:
formula 1
\[S^2= \frac {SC_{totales}}{gl}\] Donde \[SC= Suma~de~cuadrados\\ S^2=Varianza~ total\\ gl= Grados~de~libertad~ (n-1)\]
y tambien sabemos que:
Formula 2
\[SC_{totales}= SC_{trt}~ + ~ SC_{error}\] Entonces reemplazamos la formula dos en la formula uno obteniendo este resultado
formula 3
\[S^2= \frac {SC_{trt}~ + ~ SC_{error}}{gl}\] y despejamos la Suma de Cuadrados del error
\[SC_{error}= (S^2*gl)~ - SC_{trt)}\]
SCerror=(823*71)-6000;SCerror## [1] 52433\[SC_{error}=52433\] como ya tenemos las sumas de cuadrados y los grados de libertad podemos obtener las varianzas o Cuadrados Medios (CM) de los tratamientos y del error
CMtrt<- SCtrt/gl_trt;CMtrt## [1] 750CMerror<- SCerror/gl_error; CMerror## [1] 832.2698\[CM_{trt}=750\\ CM_{error}=832.26\]
ya con los dos cuadrados medios podemos calcular el COCIENTE F (CF) así:
CF <- CMtrt/CMerror; CF## [1] 0.90115f_cal<- CF\[Valor~ F= 0.90115\] Todos los datos obtenidos anteriormente se aprecian mejor en la siguiente tabla.
library(readxl)
punto_1 <- read_excel("punto1.xlsx")## New names:
## * `` -> ...1df<- data.frame(punto_1)
df## ...1 SC df CM F_cal
## 1 Between 6000 8 750.00 o.9011
## 2 Within 52433 63 832.26 <NA>
## 3 total 58433 71 NA <NA>¿Qué nos dice este valor de F?
lo primero que nos damos cuenta es que este valor es menor a 1, lo que indica que la variabilidad entre las repeticiones es mayor que la variabilidad entre tratamientos, NO PODEMOS CONCLUIR NADA CON RESPECTO A LOS TRATAMIENTOS YA QUE NO SON LA PRINCIPAL FUENTE DE VARIABILIDAD.
a.2
si el F tabulado es 2.8. ¿qué puede decirse acerca de la hipótesis nula de igualdad de los promedios del índice en todas las condiciones de tratamiento?
f_tab<- 2.8
x <- seq( -4, 4, by = 0.1)
y <- dnorm( x )
plot( function(x) df( x, df1 = 8, df2= 63), 0, 4, ylim = c( 0, 1 ),
col = "red", type = "l", lwd = 2,
main = "Función densida F de Fisher df1=8 y df2=63" )
abline(v=f_tab,col="blue")
text(2.9,0.2,"f_tab")
abline(v=f_cal, col = "green")
text(1.0,0.2,"f_cal")
text(3.5, 0.5, "zona de rechazo")
text(2, 0.5, "zona de no rechazo")
Como se observa en la grafica anterior nuestro F calculado cae en la zona de no rechazo de la hipótesis nula, con lo que pdrímos conclurir que los tratamientos tienenigualdad de promedios estadisticamente hablando, PERO estos datos no pueden ser tomados para desiciones agronómicas porque el valor de F fue menor de 1.
a.3
¿Vale la pena comparar las medias de tratamientos a posteriori del ANOVA (prueba de Tukey)?
NO
la prueba de tukey nos permite hacer comparaciones entre los tratamientos para ver cual fue el tratamiento causante de la diferencia, pero como en este caso los causante de la variabilidad son las repeticiones y no los tratamientos, una prueba de Tukey no es util.
Punto 2
Antes de hilar el algodón, éste debe ser procesado para eliminar las materias extrañas y la humedad. El limpiador de pelusas más común es el limpiador de pelusas tipo sierra de batería controlada. Aunque el limpiador de pelusas de motor de sierra (M1) es uno de los más efectivos, también es uno de los limpiadores que causa más daño a las fibras de algodón. Un investigador del algodón diseñó un estudio para comparar cuatro alternativas de limpieza de las fibras de algodón: M2, M3, M4 y M5. Los métodos M2 y M3 son mecánicos, mientras que los métodos M4 y M5 son una combinación mecánica y química. El investigador quiso tener en cuenta el impacto de los diferentes cultivadores en el proceso y, por lo tanto, obtuvo fardos de algodón de seis diferentes granjas algodoneras. Las granjas fueron consideradas como bloques en el estudio. Después de una limpieza preliminar de algodón, los seis fardos fueron mezclados a fondo, y luego fue procesada una igual cantidad de algodón por cada uno de los cinco métodos de limpieza de pelusas. Las pérdidas en peso (en kg) después de la limpieza las fibras de algodón se dan en la siguiente tabla. Durante el procesamiento de las muestras de algodón, las mediciones de la granja 1 procesada por el limpiador M1 se perdieron.
Tabla de granjero vs metodo
| Método | 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|---|
| M1 | * | 6.75 | 13.05 | 10.26 | 8.01 | 8.42 |
| M2 | 5.54 | 3.53 | 11.20 | 7.21 | 3.24 | 6.45 |
| M3 | 7.67 | 4.15 | 9.79 | 8.27 | 6.75 | 5.50 |
| M4 | 7.89 | 1.97 | 8.97 | 6.12 | 4.22 | 7.84 |
| M5 | 9.27 | 4.39 | 13.44 | 9.13 | 9.20 | 7.13 |
| Mean | 7.593 | 4.158 | 11.290 | 8.198 | 6.280 | 7.068 |
- Realice el ANOVA para este diseño recordando que es un caso desbalanceado. Concluya sobre el resultado de la tabla del ANOVA obtenida. (¿Afecta el orden de colocación de los efectos del modelo dentro del software R? Verifique si la tabla del ANOVA cambia): en R: lm()
library(readxl)
#se crea un dataframe y se llama el excel con los datos
dalgodon <- read_excel("tablae2.xlsx")
# View(tablae2)
# Se crean factores para realizar el anova
dalgodon$GRANJERO = as.factor(dalgodon$GRANJERO)
dalgodon$METODO = as.factor(dalgodon$METODO)
library(daewr)## Warning: package 'daewr' was built under R version 4.0.3## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.designMODELO 1
# se utiliza la funcion lm
mod1 <- lm( PESO ~ GRANJERO * METODO, data = dalgodon )
anova( mod1 )## Warning in anova.lm(mod1): ANOVA F-tests on an essentially perfect fit are
## unreliable## Analysis of Variance Table
##
## Response: PESO
## Df Sum Sq Mean Sq F value Pr(>F)
## GRANJERO 5 138.30 27.6608
## METODO 4 49.12 12.2799
## GRANJERO:METODO 19 26.23 1.3805
## Residuals 0 0.00table (dalgodon$GRANJERO, dalgodon$METODO)##
## M1 M2 M3 M4 M5
## 1 1 1 1 1 1
## 2 1 1 1 1 1
## 3 1 1 1 1 1
## 4 1 1 1 1 1
## 5 1 1 1 1 1
## 6 1 1 1 1 1No se puede utilizar este modelo por que solo hay una respuesta por bloque y por un factor como se observa en la tabla
MODELO 2
mod2 <- lm( PESO ~ GRANJERO + METODO, data = dalgodon )
anova( mod2 )## Analysis of Variance Table
##
## Response: PESO
## Df Sum Sq Mean Sq F value Pr(>F)
## GRANJERO 5 138.30 27.6608 20.0365 5.57e-07 ***
## METODO 4 49.12 12.2799 8.8951 0.0003186 ***
## Residuals 19 26.23 1.3805
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Este metodo es el correcto ya que primero se coloca el bloque y luego la razon experimental entonces se bloquea primero el Granjero por que la variable de interes es el metodo
MODELO 3 En este modelo se intercalaron las variables metodo y granero
mod3 <- lm( PESO ~ METODO + GRANJERO, data = dalgodon )
anova( mod3 )## Analysis of Variance Table
##
## Response: PESO
## Df Sum Sq Mean Sq F value Pr(>F)
## METODO 4 47.763 11.9407 8.6494 0.0003754 ***
## GRANJERO 5 139.661 27.9322 20.2331 5.163e-07 ***
## Residuals 19 26.230 1.3805
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se encontraron diferencias significativas en los modelos 2 y 3
- Estimar el valor de la observación usando el promedio de los datos para los cinco granjeros del mismo método M1 y luego realice el análisis de varianza para probar las diferencias en las pérdidas medias de peso para los cinco métodos de limpiado de las fibras de algodón. Compare este resultado con el caso desbalanceado (de ser posible).
# Crea el data frame nuevamente para hacer uso de los datos
dalgodon2 <- dalgodonSe reemplaza el dato G1 para el metodo 1, en su lugar se coloca el dato de la media
dalgodon2$PESO[1] = 7.593MODELO 4
mod4<- lm( PESO ~ GRANJERO + METODO, data = dalgodon2 )
anova( mod4 )## Analysis of Variance Table
##
## Response: PESO
## Df Sum Sq Mean Sq F value Pr(>F)
## GRANJERO 5 138.331 27.6662 18.4550 7.003e-07 ***
## METODO 4 45.367 11.3418 7.5656 0.0006986 ***
## Residuals 20 29.982 1.4991
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Al comparar el modelo 2 que se dijo era el indicado y el modelo 4 se observan diferencias por lo que no utilizaria este modelo.
- Probemos nuestro ingenio y capacidad de trabajo en R: Proponga un promedio basado en la medida de M1 con los datos disponibles, pero sumando una cantidad delta de modo que el coeficiente de variación entre las seis mediciones sea menor al 20%. Es decir,
Punto 3
Use la función de R para generar de la distribución uniforme unos datos de carbono orgánico del suelo medida a 5 cm y 10 cm de profundidad. Suponga que la medida de la capa superior osciló entre 3.0 y 3.U+0.1 y de la capa inferior osciló entre 2 y 2.T+0.2. Use expand.grid para generar una ventana de observación de 0 a 100 m para la longitud y de 0 a 200 m para la latitud. Genere 50 datos en cada capa. Use la función sort.int de R para ordenar los datos de cada capa con la opción partial=25+U dentro de la propia función sort.int. Una vez cree los datos realice algún diagrama de color (preferiblemente 3D) que permita visualizar las medidas de carbono en cada capa generadas por computadora. Compare si se encuentran diferencias en la media de carbono entre capas utilizando un nivel de confianza del 95%.
Desarrollo
Se fija la semilla para obtener los mismos datos y en este ejercicio se remplaza la U por 0 y T por 2
Se generan 50 datos uniformes con max y min para una produndidad de 5 y de 10, y se ordenan con la funcion sort.int
set.seed(2016)
co5<-runif(50,min =3.0, max = 3.1)
co5 =sort.int(co5,partial = 25)
co10 <- runif(50, min = 2, max = 2.4)
co10 =sort.int(co10,partial = 25)se crea la variable ventana y se definen los puntos para generar los 50 datos
ventana <- expand.grid(longitud= seq(0,100,25), latitud = seq(0,200,length.out = 10))Se crea el dataframe para realizar la grafica 3D
dfco <- data.frame(longitud = rep(ventana$longitud,2),
latitud = rep(ventana$latitud,2),
profundidad = rep(c(5,10),each = 50),
co = c(co5,co10))Grafica 3D
library(plotly)## Warning: package 'plotly' was built under R version 4.0.3## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutplot_ly (x=dfco$longitud, y=dfco$latitud, z=dfco$profundidad, type="scatter3d", mode="markers", color=dfco$co)## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.prueba de shapiro.test para saber la normalidad de los datos
x.testco5 <- shapiro.test(co5)
print(x.testco5)##
## Shapiro-Wilk normality test
##
## data: co5
## W = 0.94037, p-value = 0.01394x.testco10 <- shapiro.test(co10)
print(x.testco10)##
## Shapiro-Wilk normality test
##
## data: co10
## W = 0.91431, p-value = 0.001477Como se rechaza la hipotesis nula se dice que se debe realizar la prueba de wilcoxon en lugar de la T-student
var.test(co5, co10)##
## F test to compare two variances
##
## data: co5 and co10
## F = 0.049713, num df = 49, denom df = 49, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.02821105 0.08760401
## sample estimates:
## ratio of variances
## 0.04971319# no tiene varianzas iguales
wilcox.test(x = co5, y = co10, paired = TRUE, alternative = "t", conf.int = 0.95)##
## Wilcoxon signed rank test with continuity correction
##
## data: co5 and co10
## V = 1275, p-value = 7.79e-10
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
## 0.8320447 0.9004684
## sample estimates:
## (pseudo)median
## 0.8654967De acuerdo al p valor de la prueba wilcoxon se rechaza la hipotesis nula por lo que son diferentes en la media de carbono entre capas utilizando un nivel de confianza del 95%
Punto 4
- El siguiente diseño se corresponde con un factorial completo (32) en arreglo completamente al azar. Los factores y la respuesta fueron creados con el código:
Punto 5
Existe un tipo de diseño anidado (factorial incompleta) conocido como anidado escalonado (staggered nested design) y ocurre tal como se muestra en la imagen, donde se tienen dos fincas sembradas con variedades de papa solo que la finca A permite que se desarrollen las dos variedades mientras que la altitud de la finca B solo permite el desarrollo de una de ellas. Además, se tienen dos parcelas con la variedad 1 en la primera finca y solo una en el resto de las fincas.
Use la librería lme4 tal como aparece en el código abajo. La etiqueta “ue” hace referencia a la unidad experimental (parcela) utilizada, por lo que se necesita crear una columna que identifique la parcela, una que identifique la finca, otra para la variedad y otra para lo que aquí se llama test pero que hace referencia en este caso a los cuadrados de 1.5m*1.5m usados para tomar las muestras de plantas dentro de las parcelas.
Estos diseños son usados para estimar la varianza atribuible a las parcelas, a las parcelas anidadas en las fincas, y a la variedad dentro de la finca. El código presentado puede ayudar a la estimación de estas varianzas.
Use los datos que se muestran para estimar las varianzas antes descritas. Una ayuda para la solución de este problema puede encontrarse en el libro: Design and Analysis of Experiments with R de John Lawson.
library(readxl)
punto5 <- read_excel("punto5.xlsx")
# View(punto5)
df = data.frame(punto5)
df## finca variedad test parcela respuesta
## 1 1 1 1 1 9.76
## 2 1 1 1 2 10.65
## 3 1 1 1 3 6.50
## 4 1 1 1 4 8.08
## 5 1 1 1 5 7.84
## 6 1 1 1 6 9.00
## 7 1 1 1 7 12.81
## 8 1 1 1 8 10.62
## 9 1 1 1 9 4.88
## 10 1 1 1 10 9.38
## 11 1 1 1 11 5.91
## 12 1 1 1 12 7.19
## 13 1 1 1 13 7.93
## 14 1 1 1 14 3.70
## 15 1 1 1 15 4.64
## 16 1 1 1 16 5.94
## 17 1 1 1 17 9.50
## 18 1 1 1 18 10.93
## 19 1 1 1 19 11.95
## 20 1 1 1 20 4.34
## 21 1 1 2 1 9.24
## 22 1 1 2 2 7.77
## 23 1 1 2 3 6.26
## 24 1 1 2 4 5.28
## 25 1 1 2 5 5.91
## 26 1 1 2 6 8.38
## 27 1 1 2 7 13.58
## 28 1 1 2 8 11.71
## 29 1 1 2 9 4.96
## 30 1 1 2 10 8.02
## 31 1 1 2 11 5.79
## 32 1 1 2 12 7.22
## 33 1 1 2 13 6.48
## 34 1 1 2 14 2.86
## 35 1 1 2 15 5.70
## 36 1 1 2 16 6.28
## 37 1 1 2 17 8.00
## 38 1 1 2 18 12.15
## 39 1 1 2 19 10.58
## 40 1 1 2 20 5.45
## 41 1 2 1 1 11.91
## 42 1 2 1 2 10.00
## 43 1 2 1 3 8.02
## 44 1 2 1 4 9.15
## 45 1 2 1 5 7.43
## 46 1 2 1 6 7.01
## 47 1 2 1 7 11.13
## 48 1 2 1 8 14.07
## 49 1 2 1 9 4.08
## 50 1 2 1 10 6.73
## 51 1 2 1 11 6.59
## 52 1 2 1 12 5.77
## 53 1 2 1 13 8.12
## 54 1 2 1 14 3.95
## 55 1 2 1 15 5.96
## 56 1 2 1 16 4.18
## 57 1 2 1 17 11.25
## 58 1 2 1 18 9.51
## 59 1 2 1 19 16.79
## 60 1 2 1 20 7.51
## 61 2 1 1 1 9.02
## 62 2 1 1 2 13.69
## 63 2 1 1 3 7.95
## 64 2 1 1 4 7.46
## 65 2 1 1 5 6.11
## 66 2 1 1 6 8.58
## 67 2 1 1 7 10.00
## 68 2 1 1 8 14.56
## 69 2 1 1 9 4.76
## 70 2 1 1 10 6.99
## 71 2 1 1 11 6.55
## 72 2 1 1 12 8.33
## 73 2 1 1 13 7.43
## 74 2 1 1 14 5.92
## 75 2 1 1 15 5.88
## 76 2 1 1 16 5.24
## 77 2 1 1 17 11.14
## 78 2 1 1 18 12.71
## 79 2 1 1 19 13.08
## 80 2 1 1 20 5.21library(daewr)
mod2<-aov(respuesta ~ parcela + parcela:finca + parcela:finca:variedad, data = df)
summary(mod2)## Df Sum Sq Mean Sq F value Pr(>F)
## parcela 1 5.3 5.289 0.608 0.438
## parcela:finca 1 3.1 3.080 0.354 0.553
## parcela:finca:variedad 1 5.2 5.185 0.597 0.442
## Residuals 76 660.6 8.692library(lme4)## Warning: package 'lme4' was built under R version 4.0.3## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:daewr':
##
## cakemodr3 <- lmer( respuesta ~ 1 + (1|parcela) + (1|parcela:finca)+
+ (1|parcela:finca:variedad), data = df)## boundary (singular) fit: see ?isSingularsummary(modr3)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## respuesta ~ 1 + (1 | parcela) + (1 | parcela:finca) + +(1 | parcela:finca:variedad)
## Data: df
##
## REML criterion at convergence: 326
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92792 -0.39924 0.00919 0.43823 1.65354
##
## Random effects:
## Groups Name Variance Std.Dev.
## parcela:finca:variedad (Intercept) 1.2309 1.1094
## parcela:finca (Intercept) 0.0000 0.0000
## parcela (Intercept) 7.0122 2.6481
## Residual 0.8789 0.9375
## Number of obs: 80, groups:
## parcela:finca:variedad, 60; parcela:finca, 40; parcela, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8.2368 0.6188 13.31
## convergence code: 0
## boundary (singular) fit: see ?isSingularCon los datos obtenidos se saca la varianza total, esta es igual a 100 factor de la varianza de parcela, dividido entre la suma de varianza de parcela, varianza de parcela finca variedad y los residuales.
\[ (100*7.0122)/(7.0122+1.2309+0.8789)) = 77 \% \] De estos resultados, vemos que el 77% de la variación total se debe a la variabilidad entre parcelas, mientras que dentro de las parcelas (parcelas:finca) la variabilidad de caja a caja es insignificante.
Punto 6
Se tienen unos datod de potasio de muestras de suelos medidas en 8 diferentes laboratorios, Compare descriptivamente los datos.
- los datos son los siguientes:
library(asbio)## Warning: package 'asbio' was built under R version 4.0.3## Loading required package: tcltkdata(K)
soil<- K
soil## K lab
## 1 296 B
## 2 260 B
## 3 341 B
## 4 359 B
## 5 323 B
## 6 321 B
## 7 287 B
## 8 413 B
## 9 335 B
## 10 315 D
## 11 330 D
## 12 326 D
## 13 354 D
## 14 266 D
## 15 348 D
## 16 343 D
## 17 284 D
## 18 324 D
## 19 351 E
## 20 302 E
## 21 395 E
## 22 357 E
## 23 400 E
## 24 187 E
## 25 376 E
## 26 283 E
## 27 198 E
## 28 327 F
## 29 354 F
## 30 308 F
## 31 274 F
## 32 324 F
## 33 305 F
## 34 347 F
## 35 297 F
## 36 305 F
## 37 326 G
## 38 301 G
## 39 316 G
## 40 312 G
## 41 297 G
## 42 280 G
## 43 300 G
## 44 319 G
## 45 286 G
## 46 218 H
## 47 280 H
## 48 241 H
## 49 226 H
## 50 243 H
## 51 199 H
## 52 205 H
## 53 225 H
## 54 227 H
## 55 338 I
## 56 303 I
## 57 341 I
## 58 311 I
## 59 355 I
## 60 269 I
## 61 284 I
## 62 279 I
## 63 339 I
## 64 359 J
## 65 318 J
## 66 313 J
## 67 352 J
## 68 334 J
## 69 356 J
## 70 342 J
## 71 299 J
## 72 353 J¿qué sabemos?
Partimos suponiendo que las muestras son indepedientes, osea que la manera de recolección fue totalmente aleatoria, pero cuando estemos comprobando los supuestos también verificaremos esta hipótesis. Se sabe que las misma muestra se distribuyó a los 8 laboratorios disponibles, por eso sabemos que son muestras pareadas. Tambien al ser las misma muestras y como se espera que todos los laboratorios manejen los mismos metodos para medir el potasio del suelo, esperamos los resultados no varíen demasiado cuando comparamos un laboratorio con otro.
primero vamos a obtener un analisis descriptivo por cada uno de los lboratorios
library(psych)##
## Attaching package: 'psych'## The following object is masked from 'package:asbio':
##
## skew## The following objects are masked from 'package:ggplot2':
##
## %+%, alphadescribe(soil, ranges=T)## vars n mean sd median trimmed mad min max range skew kurtosis
## K 1 72 307.79 48.40 314.0 311.1 43.74 187 413 226 -0.54 0.04
## lab* 2 72 4.50 2.31 4.5 4.5 2.97 1 8 7 0.00 -1.29
## se
## K 5.70
## lab* 0.27describeBy(soil, soil$lab)##
## Descriptive statistics by group
## group: B
## vars n mean sd median trimmed mad min max range skew kurtosis se
## K 1 9 326.11 44.41 323 326.11 40.03 260 413 153 0.41 -0.7 14.8
## lab* 2 9 1.00 0.00 1 1.00 0.00 1 1 0 NaN NaN 0.0
## ------------------------------------------------------------
## group: D
## vars n mean sd median trimmed mad min max range skew kurtosis se
## K 1 9 321.11 29.26 326 321.11 25.2 266 354 88 -0.68 -1.02 9.75
## lab* 2 9 2.00 0.00 2 2.00 0.0 2 2 0 NaN NaN 0.00
## ------------------------------------------------------------
## group: E
## vars n mean sd median trimmed mad min max range skew kurtosis
## K 1 9 316.56 80.35 351 316.56 72.65 187 400 213 -0.54 -1.44
## lab* 2 9 3.00 0.00 3 3.00 0.00 3 3 0 NaN NaN
## se
## K 26.78
## lab* 0.00
## ------------------------------------------------------------
## group: F
## vars n mean sd median trimmed mad min max range skew kurtosis se
## K 1 9 315.67 25.05 308 315.67 23.72 274 354 80 0.05 -1.22 8.35
## lab* 2 9 4.00 0.00 4 4.00 0.00 4 4 0 NaN NaN 0.00
## ------------------------------------------------------------
## group: G
## vars n mean sd median trimmed mad min max range skew kurtosis se
## K 1 9 304.11 15.37 301 304.11 22.24 280 326 46 -0.14 -1.51 5.12
## lab* 2 9 5.00 0.00 5 5.00 0.00 5 5 0 NaN NaN 0.00
## ------------------------------------------------------------
## group: H
## vars n mean sd median trimmed mad min max range skew kurtosis se
## K 1 9 229.33 23.89 226 229.33 22.24 199 280 81 0.74 -0.32 7.96
## lab* 2 9 6.00 0.00 6 6.00 0.00 6 6 0 NaN NaN 0.00
## ------------------------------------------------------------
## group: I
## vars n mean sd median trimmed mad min max range skew kurtosis se
## K 1 9 313.22 31.4 311 313.22 41.51 269 355 86 -0.09 -1.81 10.47
## lab* 2 9 7.00 0.0 7 7.00 0.00 7 7 0 NaN NaN 0.00
## ------------------------------------------------------------
## group: J
## vars n mean sd median trimmed mad min max range skew kurtosis se
## K 1 9 336.22 21.61 342 336.22 20.76 299 359 60 -0.46 -1.53 7.2
## lab* 2 9 8.00 0.00 8 8.00 0.00 8 8 0 NaN NaN 0.0para una mejor interpretación se observarán solo los valores de media, mediana y desviación estándar.
mean<-tapply(soil$K, soil$lab,mean);mean## B D E F G H I J
## 326.1111 321.1111 316.5556 315.6667 304.1111 229.3333 313.2222 336.2222sd<-tapply(soil$K, soil$lab,sd);sd## B D E F G H I J
## 44.40564 29.25510 80.35097 25.04995 15.37404 23.89037 31.39577 21.60890median<- tapply(soil$K, soil$lab,median);
u_u<-data.frame(mean,sd,median)
library(DT)## Warning: package 'DT' was built under R version 4.0.3datatable(u_u, class = 'cell-border stripe',filter = 'top', options = list(
pageLength = 5, autoWidth = TRUE))ahora, se va a relaizar nuestro primer gráfico que nos va a permitir analizar de manera más sencilla nuestros datos.
library(ggplot2)
ggplot(soil, aes(x = lab, y = K))+
geom_boxplot()
ahora que tenemos una apreciación visual podemos ver que: - los las muestras entre laboratorios no son similares - la distribución de los datos en los boxplot no parecen tener una distribución normal - la varianza de los datos entre laboratorios tambien varían demasiado
¿Qué prueba usar para evaluar si las medias son iguales o no?
ya que observamos que las medias son distintas podemos evaluar a través de la prueba de Kruskal- Wallis si la variable respuesta es la misma en todas las poblaciones valoradas (los ocho laboratorios), osea si las variables pertenen a una misma distribución, esto a través de rangos.
\[H_o= la~variable~respuesta~es ~la~misma~para~todos~los~laboratorios \\ H_a= Al~menos~un~laboratorio~presenta~una~variable~distinta\]
kr<- kruskal.test(soil$K,soil$lab)
kr$statistic## Kruskal-Wallis chi-squared
## 24.48198a partir de qué valor podemos recharaz o aprobar la Ho?
qchisq(0.05, 7, lower.tail = F)## [1] 14.06714dado que nuestro valor de K mayor a nuestro valor estimado, rechazamos nuestra H_O, y concluimos que todos nuestros datos no pertenecen a un mismo tipo de población, es decir, las medias de los rangos entre laboratorios no son las mismas.
Ahora vamos a comparar las diferencias pareadas para ver dónde están las diferencias así:
PMCMR::posthoc.kruskal.nemenyi.test(soil$K~soil$lab)## Warning in posthoc.kruskal.nemenyi.test.default(c(296, 260, 341, 359, 323, :
## Ties are present, p-values are not corrected.##
## Pairwise comparisons using Tukey and Kramer (Nemenyi) test
## with Tukey-Dist approximation for independent samples
##
## data: soil$K by soil$lab
##
## B D E F G H I
## D 1.0000 - - - - - -
## E 1.0000 1.0000 - - - - -
## F 0.9999 0.9999 0.9998 - - - -
## G 0.9324 0.9324 0.9222 0.9943 - - -
## H 0.0098 0.0098 0.0087 0.0397 0.2764 - -
## I 0.9993 0.9993 0.9989 1.0000 0.9984 0.0600 -
## J 0.9893 0.9893 0.9916 0.9051 0.4405 0.0003 0.8461
##
## P value adjustment method: none¿qué podemos concluir?
El laboratorio que sus media de rangos difiere con las medias de rangos de la mayoría de los otros laboratorios (6 de 7) es el laboratorio H, si comparamos los resultados se los otros laboratior no muestran diferencias significativas estadísticamente hablando.