PARCIAL DISENIO DE EXPERIMENTOS

Universidad Nacional de Colombia/Facultad de Ciencias Agrarias/Dpto. de Agronomía Diseño de experimentos/Parcial I

Apellidos y Nombres: Pablo Andrés Cárdenas Rodríguez, Eduardo Ducuara Alape Cédula: 1058038901, 1010244260 Dirección de RPubs para la revisión:

Cédula usada para la creación de datos (digitos faltantes U=1,T=5)= 1058038901

a- En un estudio conducido en ambiente controlado se tuvieron 72 macetas, cada una con una planta a la que a cierta edad se le midió el contenido de clorofila (índice de clorofila) con un sensor (SPAD). El total de macetas se correspondió con 9 tratamientos asociados a estrés hídrico. Se sabe que la varianza de las 72 observaciones es 813. Con esta información complete la tabla del ANOVA que se muestra a continuación:

	Suma de cuadrados	gl	Media cuadrada	\(F\)
Between	6000	-	-	-
Within	-	-	-
Total	-	-

Solución

Datos suministrados

\[n=72\\Tratamientos(trt)=9\\Repeticiones(rep)=8\\Varianza=813\\Suma~Cuadrados_{trt}=6000\\Grados~ de~libertad~tratamientos(gl_{trt})=trt-1=8\\Grados~de~libertad~error(gl_{error})=trt*(rep-1)=9(7)=63\\Grados~de~libertad~totales(gl_{total})=n-1=71\]

Para completar la tabla ANOVA es necesario deducir la suma de cuadrados de las repeticiones \(SC_{rep}\), y/o la suma de cuadrados del error\(SC_{error}\). para ello podemos usar las fórmula propuestas por R. Fisher:

-Fórmula (1)

\[S^2=\frac{SC_{total}}{gl}\] \[*Sean:\\SC_{total}=Suma~de~cuadrados~Total\\S^2=Varianza~Total\\gl=Grados~de~libertad~Totales\]

Para la fórmula 1 es posible hallar la suma de cuadrados totales usando:

Fórmula (2)

\[SC_{total}=SC_{trt}+SC_{error}\] \[*Donde:\\SC_{total}=Suma~de~cuadrados~Total\\SC_{trt}=Suma~de~cuadrados~Tratamientos\\SC_{error}=Suma~de~cuadrados~Error\]

Y reemplazando la fórmula (2) en (1):

Fórmula (3)

\[S^2=\frac{SC_{trt}+SC_{error}}{gl_{tot}}\]

Y despejando \(SC_{error}\) obtenemos

\[SC_{error}=(S^2*gl_{tot})-SC_{trt}\]

SCerror=(813*71)-6000
SCerror

## [1] 51723

\[SC_{error}=51723\] Al tener la suma de cuadrados y los grados de libertad es posible obtener las varianzas o Cuadrados medios de tratamientos \(CM_{trt}\) y del error \(CM_{error}\)

S2=823
SCtrt=6000
gl_trt=8
gl_error=63
gl_totales= 71

CMtrt=SCtrt/gl_trt
CMtrt

## [1] 750

CMerror=SCerror/gl_error
CMerror

## [1] 821

\[CM_{trt}=750\\CM_{error}=821\] Para hallar el cociente \(F\) usamos: \[F=\frac{CM_{trt}}{CM_{error}}\]

F=CMtrt/CMerror
F

## [1] 0.9135201

\[F~valor=0.9135201\]

Tabla ANOVA

	Suma de cuadrados	gl	Media cuadrada	\(F\)
Entre trt	6000	8	750	0.9135201
Entre Rep	51723	63	821
Total	57723	71

Conclusión Tabla ANOVA

El F valor calculado es menor a 1, por lo que la variabilidad entre las repeticiones es mayor que la debida a los tratamientos, habría que evaluar que tan diferentes son comparando el estadistico de prueba F calculado.

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 (use el p valor así como el cociente F calculado de la tabla para concluir)?

Hipotesis nula

\[H_0: \mu_{trt1}=\mu_{trt2}=\mu_{trt3}=\mu_{trt4}=\mu_{trt5}=\mu_{trt5}=\mu_{trt6}=\mu_{trt7}=\mu_{trt8}=\mu_{trt9} \\H_a: \mu_{trt}~no~son~iguales ~(Al~menos~uno~diferente)\]

Usando gráfico F de Fisher con los datos de la tabla ANOVA

#Valores de F
F_tab= 2.8
F_cal= F

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 = "gold3", type = "l", lwd = 2,
      main = "Función de densidad F de Fisher df1=8 y df2=63" )
abline(v=F_tab,col="salmon1", lwd = 2)
text(2.9,0.1,"F_tab")
abline(v=F_cal, col = "aquamarine4", lwd = 2)
text(1.0,0.1,"F_cal")
text(3.5, 0.3, "Zona de rechazo")
text(0.95, 0.3, "Zona de no rechazo")
text(F_tab+0.2,0.01, "α=5%")
pf(q=F_cal, df1=8, df2=63, lower.tail=F) #Calculando la probabilidad con respecto a la distribución F usando datos hallados en la tabla Anova

## [1] 0.5113406

segments(x0=F_cal,y0=0.8, x1=4, y1=0.8, col="blue1")
text(F_cal+.62, 0.85, "p.valor=0.5113406")

-¿vale la pena comparar las medias de tratamientos a posteriori del ANOVA (prueba de Tukey)?

La hipótesis nula de que las medias de los grupos son estadísticamente iguales no se rechaza, necesitamos un valor F alto, cosa que no ocurre y es explicito en la gráfica y comparando el F calculado y el tabulado. Para rechazar se necesitaría un valor mayor a 2.8 en el coeficiente F. El p.valor calculado también es mucho mayor a 0.05 que sería un buen indicativo para no rechazar la igualdad entre tratamientos.

En este caso NO sería necesario aplicar esta prueba prueba de Tukey, dado que no se rechazó la hipótesis nula usando la tabla Anova. Aunque debería evaluarse que tan importante es la variablidad debida a las repeticiones y si es concluyente en el analisis de los datos.

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 del 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.

			Granjero
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
Media	7.593	4.158	11.290	8.198	6.280	7.068

a-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)
Algodon<- read_excel("C:/Users/Equipo/Downloads/Algodón.xlsx") #dataframe nuevo con columnas perdida(peso),método y granjero y se importa
head(Algodon)

## # A tibble: 6 x 4
##     `#` Metodo Perdida Granjero
##   <dbl> <chr>  <chr>      <dbl>
## 1     1 M1     <NA>           1
## 2     1 M1     6.75           2
## 3     1 M1     13.05          3
## 4     1 M1     10.26          4
## 5     1 M1     8.01           5
## 6     1 M1     8.42           6

#Se crean nuuevos factores para realizar el ANOVA
Algodon$Granjero=as.factor(Algodon$Granjero)
Algodon$Metodo=as.factor(Algodon$Metodo)

#library(daewr)

MODELO 1

# Se utiliza la función lm
modelo1<-lm(Perdida~Granjero*Metodo, data=Algodon)
anova(modelo1)

## Analysis of Variance Table
## 
## Response: Perdida
##                 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.00

No 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

modelo2 <- lm(Perdida~Metodo*Granjero, data = Algodon )
anova(modelo2)

## Warning in anova.lm(modelo2): ANOVA F-tests on an essentially perfect fit are
## unreliable

## Analysis of Variance Table
## 
## Response: Perdida
##                 Df  Sum Sq Mean Sq F value Pr(>F)
## Metodo           4  47.763 11.9407               
## Granjero         5 139.661 27.9322               
## Metodo:Granjero 19  26.230  1.3805               
## Residuals        0   0.000

Este 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

modelo3<-lm(Perdida~Metodo+Granjero, data=Algodon)
anova(modelo3)

## Analysis of Variance Table
## 
## Response: Perdida
##           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 ' ' 1

Se encontraron diferencias significativas en los modelos 2 y 3

b-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).

# Se renombra el dataframe para propósito de los datos de la parte b y evitar interferencias de nombres y posteriormente se le integra el dato faltante usando la media correspondiente al Método 1, Granjero 1

Algodón<-Algodon
mean=tapply(Algodón$Granjero, Algodón$Metodo, mean)

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

Algodón$Perdida[1]=7.593

MODELO 4

modelo4<- lm(Perdida~Granjero+Metodo, data = Algodon)
anova(modelo4)

## Analysis of Variance Table
## 
## Response: Perdida
##           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 ' ' 1

M2 y M4 presentan mayor variación,M1 presenta el menor CV pero también es el que ocasiona mayores perdidas de peso en el algodon (descartado). M3 podria ser elm recomendado. ordenar los factores parece no tener difer3ncia. Los modelos con interacción difieren de los sin interacción.

Al comparar el modelo 2 que se dijo era el indicado y el modelo 4 se observan diferencias por lo que no utilizaria este modelo

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.1+0.1 (3.2) y de la capa inferior osciló entre 2 y 2.5+0.2 (2.7). 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+1 (26) 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%.

datos a remplazar U=1 T=5

set.seed(2018) # Fijamos la semilla para asegurar la replicación de los datos
co5<-runif(50,min =3.0, max = 3.2) # Se generan 50 datos con distrib. uniforme para carbono a 5 cm de profundidad
co5 =sort.int(co5,partial = 26) # Se ordenan con la función sort.int y se usa un partial= 26
co10 <- runif(50, min = 2, max = 2.7)
co10 =sort.int(co10,partial = 26)

#Se procede a crear la variable ventana

ventana <- expand.grid(longitud= seq(0,100,25), latitud = seq(0,200,length.out = 10)) #  Se crea la ventana usando la función expand.grid  y se crea la secuencia de los datos para latitud y longitud

dfco <- data.frame(longitud = rep(ventana$longitud,2),
                   latitud = rep(ventana$latitud,2),
                   profundidad = rep(c(5,10),each = 50),
                   co = c(co5,co10))
# Se crea un dataframe para correr el grafico en 3D

Gráfico en 3D (Scatterplot3D)

library(plotly)

## 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':
## 
##     layout

plot_ly (x=dfco$longitud, y=dfco$latitud, z=dfco$profundidad, type="scatter3d", mode="markers", color = dfco$co,  colors = c('#BF382A', '#4b3621')) # Se crea el gráfico usando la función plot_ly

## 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.

Se realiza el test de Shapiro para esclarecer la normalidad de los datos trabajados

x.testco5 <- shapiro.test(co5)
print(x.testco5)

## 
##  Shapiro-Wilk normality test
## 
## data:  co5
## W = 0.93918, p-value = 0.0125

x.testco10 <- shapiro.test(co10)
print(x.testco10)

## 
##  Shapiro-Wilk normality test
## 
## data:  co10
## W = 0.93206, p-value = 0.006619

Según las pruebas de Shapiro, se rechaza la normalidad de las variables carbono orgánico medidas a los 5 y 10 cm.

Se realiza prueba t-pareada para comparar lasd las medias de carbono orgánico a 5 y 10 cm de profundidad

\(Prueba t\) para Dos muestras dependientes/pareadas

Hipótesis

\[H_0: \mu_{co5}=\mu_{co10}\\ H_a:\mu_{co5}\neq\mu_{co10}\]

Prueb_tco= t.test(dfco$co~dfco$profundidad, alternative = "t", paired = T, conf.level=0.95); Prueb_tco

## 
##  Paired t-test
## 
## data:  dfco$co by dfco$profundidad
## t = 33.846, df = 49, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.7327296 0.8252316
## sample estimates:
## mean of the differences 
##               0.7789806

ifelse(Prueb_tco$p.value<0.05, "Rechazo Ho","No rechazo Ho")

## [1] "Rechazo Ho"

Se rechaza la hipótesis nula de igualdad de medias para el carbono medida a 5 y 10 cm de profundidad, se establece una diferencia promedio entre las profundidas de 0.7789806. Para ilustrar el comportamiento se proponen los siguientes gráficos:

Gráficos de comportamiento de las medidas

# Se realiza un boxplot
boxplot(dfco$co~dfco$profundidad, main="Carbono organico a 5cm y a 10cm",size=0.05, col='aquamarine3') 
points(c(1,2), c(mean(co5), mean(co10)), pch =16 ,col="salmon") 
segments(1.5, mean(co5), 1.5, mean(co10), col = "orange")
text(1.5, 650, "Diferencia\n entre ambas\n medias", pos = 2)
text(1.50, mean(co5),'Media' ,font=8)
text(1.50, mean(co10),'Media' ,font=8)

Como concluye la prueba t- pareada es clara la diferencia entre los datos de carbono en función de la profundidad, con mayor variabilidad de datos para la profundidad de 10 cm y sin datos atípicos.

Gráfico de correlación y Coeficiente de correlación de Pearson

cor.test(dfco$co, dfco$profundidad, method="pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  dfco$co and dfco$profundidad
## t = -25.75, df = 98, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9547724 -0.9024289
## sample estimates:
##        cor 
## -0.9333989

plot(dfco$co~dfco$profundidad,pch=16,col="salmon", main = "Correlación de carbono orgánico a 5cm y a 10cm")
abline(dfco$co, dfco$profundidad)

Se presenta una asociación muy alta, casi perfecta de tipo negativo, y el grafico muestra una relación muy importante entre los grupos de datos de las profundidades (5 y 10) , pero son muy difderentes sus medias y su dispersión

Según el test de Shapiro , la prueba t - pareada y los gráficos realizados, se rechaza la hipotesis nula de igualdad de medias por lo que se debe realizar la prueba de wilcoxon

var.test(dfco$co, dfco$profundidad) # prueba de varianza carbono organivo vs profunidad

## 
##  F test to compare two variances
## 
## data:  dfco$co and dfco$profundidad
## F = 0.02786, num df = 99, denom df = 99, p-value < 2.2e-16
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.01874527 0.04140625
## sample estimates:
## ratio of variances 
##         0.02785985

var.test(dfco$co~dfco$profundidad) # prueba de varianza carbono en función de profundidad

## 
##  F test to compare two variances
## 
## data:  dfco$co by dfco$profundidad
## F = 0.084188, num df = 49, denom df = 49, p-value = 6.661e-15
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.04777465 0.14835503
## sample estimates:
## ratio of variances 
##         0.08418794

Las dos varianzas son diferentes sin importar los grados de libertad de los dos modelos usados.

wilcox.test(dfco$co, dfco$profundidad, paired = TRUE, alternative = "t", conf.int = 0.95)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  dfco$co and dfco$profundidad
## V = 0, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -4.872059 -4.719065
## sample estimates:
## (pseudo)median 
##      -4.812727

wilcox.test(dfco$co~dfco$profundidad, paired = TRUE, alternative = "t", conf.int = 0.95)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  dfco$co by dfco$profundidad
## V = 1275, p-value = 7.79e-10
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  0.7346542 0.8347585
## sample estimates:
## (pseudo)median 
##      0.7897409

Finalmente se concluye que las medias de carbono entre las capas son diferentes cuando se utiliza el nivel del 95% de confianza, de acuerdo a lo obtenido del p valor en las distintas pruebas y los gráficos propuestos.

PUNTO 4

Diseño factorial completo (3^2) en arreglo completamente al azar Este diseño está conformado por 2 factores (Dosis del insecticida y Número de aplicaciones), cada uno con 3 niveles, es decir, presenta 9 combinaciones de tratamientos cada uno se repite dos veces para un total de 18 observaciones. Los factores y la respuesta fueron creados con el código:

D <- expand.grid(Dosis_insecticida=c(3.25, 3.75, 4.25), Numero_aplicaciones=c(4,5,6), Repeticiones=paste0("r",1:6)) #Crea el diseño 3^2
D <- rbind(D, D) # Crea la estructura para dos repeticiones por tratamiento
set.seed(2020)
D <- D[order(sample(1:18)),] # Aleatoriza la estructura
class(D)

## [1] "data.frame"

D$biomasa=sort.int(rnorm(18,3,0.3),partial = 9) #Crea la respuesta 
head(D)

##    Dosis_insecticida Numero_aplicaciones Repeticiones  biomasa
## 2               3.75                   4           r1 2.708826
## 10              3.25                   4           r2 2.772692
## 16              3.25                   6           r2 2.143359
## 4               3.25                   5           r1 2.560519
## 13              3.25                   5           r2 2.708666
## 6               4.25                   5           r1 2.773705

Aqui el paso a seguir fue la creación de un data.frame con los datos de la conversión de las siguientes variables, (variable a factor) y (respuesta a vector):

Dosis_insecticida1=as.factor(D$Dosis_insecticida)
Numero_aplicaciones1 = as.factor(D$Numero_aplicaciones)
Biomasa1 = as.vector(D$biomasa)
D1<-data.frame(Dosis_insecticida1,Numero_aplicaciones1,Biomasa1)
str(D1)

## 'data.frame':    18 obs. of  3 variables:
##  $ Dosis_insecticida1  : Factor w/ 3 levels "3.25","3.75",..: 2 1 1 1 1 3 2 2 2 2 ...
##  $ Numero_aplicaciones1: Factor w/ 3 levels "4","5","6": 1 1 3 2 2 2 3 3 2 2 ...
##  $ Biomasa1            : num  2.71 2.77 2.14 2.56 2.71 ...

Para tener una base en la que nos podamos apoyar para identificar el comportamiento realizamos un análisis descriptivo de los datos y revisar en cual se afecta menos el crecimiento del cultivo:

lattice::bwplot(D1$Biomasa1~D1$Numero_aplicaciones1|D1$Dosis_insecticida1)

Del gráfico se puede interpretar que las plantas generan mayor biomasa es cuando tinenen la dosis 4.25 con 4 aplicaciones por lo tanto hasta elmomento se podría considerar el tratamiento 4.25 con 4 aplicaciones como una buena opción para el manejo de malezas sin tener mayor afectación en el cultivo.

Para visualizar el diseño:

library(collapsibleTree)

collapsibleTree(D,hierarchy=c("Dosis_insecticida", "Numero_aplicaciones","Repeticiones"))

a: Escriba el Modelo del diseño (completamente especificado)

\[y_{ijk}=μ+α_i+β_j+(αβ)_{ij}+ϵ_{ijk}\]

\[i:1⋯3\] \[j:1,3\] \[k:1⋯6\]

\(Donde:\)

\(y_{ijk}\): Es la ikj−esima biomasa para el i−esimo nivel del factor 1 (Dosis_insecticida) y el j−esimo nivel del factor 2 (Numero_aplicaciones).

μ: Es la media general de Biomasa

\(α_i\): Es el efecto del i−esimo nivel del factor 1 (Dosis_insecticida).

\(β_j\): Es el efecto del j−esimo nivel del factor 2 (Numero_aplicaciones).

\((αβ)_{ij}\): Es la interacción del i−esimo nivel del factor 1 (Dosis_insecticida) con el j−esimo nivel del factor 2 (Numero_aplicaciones).

\(ϵ_{ijk}\): Es el error residual para el i−esimo nivel del factor 1 (Dosis_insecticida), j−esimo nivel del factor 2 (Numero_aplicaciones) y la k−esima repetición.

Plantear Hipótesis

El efecto de la interacción entre los dos factores es nulo para todo i,j. \[H_{01}:(αβ)_{ij}=0 ∀i.j\]
El efecto del F1, en este caso el efcto de la la Dosis del insecticida es nulo. \(H_{02}:α_i=0\) \(i=1⋯3\)
El efecto del F2, en este caso el efecto del número de aplicaciones es nulo. \(H_{03}:β_j=0\) \(j=1,3\)

Realice el Anova para este diseño y de ser necesario realice la prueba de comparaciones de medias para los efectos principales de F1: dosis de un insecticida que se sospecha tiene un efecto de disminución del crecimiento (biomasa) y F2: número de aplicaciones durante el desarrollo del cultivo.

ANOVA

modelo1 <- aov(D1$Biomasa1~(D1$Dosis_insecticida1*D1$Numero_aplicaciones1))
summary(modelo1)

##                                               Df Sum Sq Mean Sq F value Pr(>F)
## D1$Dosis_insecticida1                          2 0.4161 0.20804   1.979  0.194
## D1$Numero_aplicaciones1                        2 0.3426 0.17128   1.630  0.249
## D1$Dosis_insecticida1:D1$Numero_aplicaciones1  4 0.1635 0.04087   0.389  0.812
## Residuals                                      9 0.9459 0.10510

De acuerdo con el análisis de varianza no existe interacción entre las dosis del insecticida y el numero de aplicaciones, observando el F-value, este presenta valores mayores que 1 en dos factores por lo que se supondría que la principal fuente de variación son los tratamientos y no las replicaciones, lo que sería ideal en este caso.

En cuanto a la segunda y tercera hipótesis, se dice que el efecto de los dos factores (Dosis y número de aplicaciones) es nulo sobre la biomasa, en este caso se observan que los p_valores son mayores que 0.05 lo cual nos da a decir que no es tan significativo el resultado por ende no se van a rechazar las hipotesis.

Prueba de comparación de medias: Se hizo esta prueba por no haber interacción segun los datos arrojados por el ANOVA.

modelo1b <- aov(D1$Biomasa1~D1$Dosis_insecticida1+D1$Numero_aplicaciones1)
TukeyHSD(modelo1b)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = D1$Biomasa1 ~ D1$Dosis_insecticida1 + D1$Numero_aplicaciones1)
## 
## $`D1$Dosis_insecticida1`
##                 diff         lwr       upr     p adj
## 3.75-3.25 0.27263391 -0.17270866 0.7179765 0.2738312
## 4.25-3.25 0.35602957 -0.08931299 0.8013721 0.1261203
## 4.25-3.75 0.08339567 -0.36194690 0.5287382 0.8751713
## 
## $`D1$Numero_aplicaciones1`
##            diff        lwr       upr     p adj
## 5-4 -0.24371366 -0.6890562 0.2016289 0.3479341
## 6-4 -0.32456967 -0.7699122 0.1207729 0.1711335
## 6-5 -0.08085601 -0.5261986 0.3644866 0.8821270

Haciendo la comparación de las medias del efecto de las dosis del insecticida sobre la biomasa del cultivo se observa que el p-valor ajustado en todos los casos es mayor a 0.05, por lo que podemos decir que se tendría el mismo efecto al usar diferentes dosis. Hasta este momento se va contradiciendo el análisis descriptivo hecho anteriormente, en el que con la dosis 4.25 las plantas presentan mayor biomasa, y también las comparaciones de medias del efecto del número de aplicaciones sobre la biomasa, estas tienen el mismo comportamiento que las dosis del insecticida y por tanto, su efecto sobre la biomasa es igual.

Revisión de supuestos * Prueba de normalidad de los residuales \[H_O=Residuales~normales\] \[H_O=Residuales ~no~ normales\]

#Prueba de normalidad de los residuales
residuales=modelo1b$residuals; residuales

##           1           2           3           4           5           6 
## -0.48411735 -0.14761743 -0.45238058 -0.11607671  0.03207085 -0.25891959 
##           7           8           9          10          11          12 
## -0.09802343 -0.03590285 -0.05094965  0.41038930  0.10233013  0.23758720 
##          13          14          15          16          17          18 
##  0.21133006  0.25860398 -0.01651421  0.44641666 -0.07578647  0.03756007

shres=shapiro.test(residuales)
ifelse(shres$p.value<0.05, "Rechazo Ho","No rechazo Ho")

## [1] "No rechazo Ho"

Prueba de homocedasticidad de residuales

\(H_0= Varianzas ~iguales\)

\(H_0=Almenos~una~varianza~diferente\)

# Se realiza la prueba de homocedasticidad de residuales  
HR=bartlett.test(residuales~D1$Dosis_insecticida1)
ifelse(HR$p.value<0.05, "Rechazo Ho","No rechazo Ho")

## [1] "No rechazo Ho"

Supuesto de independencia de residualeS

# Grafico de independencia de residuales
plot(residuales, pch=16, col="red",  main="Gráfico de residuales")

Como se puede ver en el gráfico no se evidencia alguna tendencia que tenga una función.

Por lo que visualizaremos la tabla de medias, para comparar nuevamente el efecto de los tratamientos sobre el crecimiento del cultivo (biomasa):

t <- tapply(Biomasa1, list(Dosis_insecticida1, Numero_aplicaciones1), mean)
## addmargins(t, FUN= mean)

De lo anterior se confirma que el mejor tratamiento sigue siendo el de la dosis 4,25 con 4 aplicaciones ya que afecta menos el crecimiento del cultivo (mayor biomasa) por lo que sigue sin existir interacción.

Gráfico de interacción

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
D1 %>%
  group_by(Dosis_insecticida1,Numero_aplicaciones1) %>%
  summarise(medias_biomasa = mean(Biomasa1)) -> mediasDQ

## `summarise()` regrouping output by 'Dosis_insecticida1' (override with `.groups` argument)

ggplot(mediasDQ)+
  aes(x=Dosis_insecticida1,y=medias_biomasa,col = Numero_aplicaciones1) +
  geom_line(aes(group = Numero_aplicaciones1))+
  geom_point()

Segun lo que acaba de pasar en la grafica anterior nos contradice lo que veniamos diciendo anteriormente de que no había interacción entre los factores Dosis del insecticida y Número de aplicaciones, en este gráfico se observa lo contrario, al parecer hay interacciones entre las aplicaciones 4 y 5 (aunque es muy poca) con las 6 aplicaciones, pero por lo que se tiene entendido en clase no es tan confiable el resultado del grafico asi que nos seguiremos basando en el ANOVA, en conclución hay interacción entre los dos factores.

El investigador quiso colocar como covariable el contenido de arcilla(expansible) en el suelo utilizado en cada unidad experimental. Genere unos datos con la distribución uniforme cuya medida oscile entre 0.20 y 0.40 , ordene estas medidas en forma decreciente y meta dentro del análisis esta variable. Especifique nuevamente el modelo y realice el análisis de covarianza respectivo ¿se justifica el uso de la covariable? Construya nuevamente el gráfico de interacción y compare con el caso sin covariable (discuta el resultado). Revise en internet los supuestos que deben tener las covariables para ser utilizadas en el modelo. ¿Se está incumpliendo en nuestros datos alguno de los supuestos necesarios? Revise los supuestos sobre los residuales tanto del ANOVA como del ANCOVA ¿qué puede percibir? ¿recomendaría el uso de arcillas para minimizar el efecto sobre el contenido de biomasa que puede ocasionar el uso del insecticida?

Covariable

set.seed(2020)
arcilla<-sort.default(runif(18, 0.2, 0.4), decreasing = TRUE) #covarable
arcilla

##  [1] 0.3921445 0.3652331 0.3528828 0.3487672 0.3307115 0.3293806 0.3240412
##  [8] 0.3237004 0.3079385 0.2953782 0.2845458 0.2818575 0.2788452 0.2786236
## [15] 0.2272194 0.2258305 0.2134769 0.2005165

ARC=data.frame(D, arcilla)
head(ARC)

##    Dosis_insecticida Numero_aplicaciones Repeticiones  biomasa   arcilla
## 2               3.75                   4           r1 2.708826 0.3921445
## 10              3.25                   4           r2 2.772692 0.3652331
## 16              3.25                   6           r2 2.143359 0.3528828
## 4               3.25                   5           r1 2.560519 0.3487672
## 13              3.25                   5           r2 2.708666 0.3307115
## 6               4.25                   5           r1 2.773705 0.3293806

Modelo del diseño incluyendo la covariable: \[y_{ijk}=μ+α_i+β_j+(αβ)_{ij}+δ(x_{ijk}−x¯+ϵ_{ijk}\] \[i:1⋯3\] \[j:1,3\] \[k:1⋯6\]

Donde: \(y_{ijk}\) Es la ikj−esima observación de biomasa para el i−esimo nivel del factor 1 (Dosis_insecticida) y el j−esimo nivel del factor 2 (Numero_aplicaciones) y la k−esima repetición.

\(μ\): Es la media general de Biomasa \(α_{i}:\) Es el efecto del i−esimo nivel del factor 1 (Dosis_insecticida).

\(β_j\): Es el efecto del j−esimo nivel del factor 2 (Numero_aplicaciones). \((αβ)_{ij}:\) Es la interacción del i−esimo nivel del factor 1 (Dosis_insecticida) con el j−esimo nivel del factor 2 (Numero_aplicaciones). \(x_{ijk}\) es la medida de la covariable que se hace para \(y_{ijk}\) \(x¯\): es la media de los valores de \(x_{ijk}\) \(δ\) :es el coeficiente de regresión que relaciona \(y_{ijk}\) con la covariable \(x_{ijk}\) \(ϵ_{ijk}\):Es el error residual para el i−esimo nivel del factor 1 (Dosis_insecticida), j−esimo nivel del factor 2 (Numero_aplicaciones) y la k−esima repetición.

Ancova Se corre el modelo con la interacción

ancova1<-aov(ARC$biomasa~ARC$arcilla+(ARC$Dosis_insecticida*ARC$Numero_aplicaciones))
summary(ancova1)

##                                               Df Sum Sq Mean Sq F value  Pr(>F)
## ARC$arcilla                                    1 0.6594  0.6594  12.238 0.00393
## ARC$Dosis_insecticida                          1 0.0681  0.0681   1.264 0.28119
## ARC$Numero_aplicaciones                        1 0.4277  0.4277   7.936 0.01454
## ARC$Dosis_insecticida:ARC$Numero_aplicaciones  1 0.0123  0.0123   0.229 0.64044
## Residuals                                     13 0.7005  0.0539                
##                                                 
## ARC$arcilla                                   **
## ARC$Dosis_insecticida                           
## ARC$Numero_aplicaciones                       * 
## ARC$Dosis_insecticida:ARC$Numero_aplicaciones   
## Residuals                                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En el análisis de covarianza se observa que la arcilla tiene un efecto significativo sobre la biomasa, esto se detecto ya que al incluir la covariable e el análisis disminuyen los grados de libertad. Es decir, se mejoró la exactitud del modelo reduciendo el error.

Ahora se corre sin la interacción (lo mejor)

ancova2<-aov(ARC$biomasa~ARC$arcilla+ARC$Dosis_insecticida+ARC$Numero_aplicaciones)
summary(ancova2)

##                         Df Sum Sq Mean Sq F value  Pr(>F)   
## ARC$arcilla              1 0.6594  0.6594  12.951 0.00291 **
## ARC$Dosis_insecticida    1 0.0681  0.0681   1.338 0.26676   
## ARC$Numero_aplicaciones  1 0.4277  0.4277   8.399 0.01168 * 
## Residuals               14 0.7128  0.0509                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Aquí el p.valor de la arcilla indica que la cantidad de arcilla si esta relacionada significativamente con la variable respuesta “Biomasa”. También se encontró significancia en el Número de aplicaciones sobre la biomasa.

Supuestos del ancova - Supuesto de linealidad entre la covariable y la biomasa.

plot(ARC$biomasa~ARC$arcilla)
Q <-lm(ARC$biomasa~ARC$arcilla,data = D)
abline(Q, col="red")

Se observa que los datos se distribuyen cerca de la linea, cumplen el supuesto.

Homogeneidad de las pendientes de regresión o interacción nula entre la covariable y la biomasa.

library(rstatix)

## 
## Attaching package: 'rstatix'

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

ARC %>% anova_test(biomasa ~ arcilla *(Dosis_insecticida+Numero_aplicaciones)+arcilla*Dosis_insecticida*Numero_aplicaciones)

## Coefficient covariances computed by hccm()

## ANOVA Table (type II tests)
## 
##                                          Effect DFn DFd      F     p p<.05
## 1                                       arcilla   1  10 11.401 0.007     *
## 2                             Dosis_insecticida   1  10  1.146 0.310      
## 3                           Numero_aplicaciones   1  10 12.090 0.006     *
## 4                     arcilla:Dosis_insecticida   1  10  6.960 0.025     *
## 5                   arcilla:Numero_aplicaciones   1  10  0.204 0.661      
## 6         Dosis_insecticida:Numero_aplicaciones   1  10  1.135 0.312      
## 7 arcilla:Dosis_insecticida:Numero_aplicaciones   1  10  0.032 0.862      
##     ges
## 1 0.533
## 2 0.103
## 3 0.547
## 4 0.410
## 5 0.020
## 6 0.102
## 7 0.003

De acuerdo a la tabla, se observa no hay interacción entre la covariable y el numero de aplicaciones pero si con la dosis ya que presenta un p-valor 0.025<0.05. No se cumple este supuesto.

Normalidad de los residuos.

G<- lm(ARC$biomasa~ARC$arcilla)
summary(G)

## 
## Call:
## lm(formula = ARC$biomasa ~ ARC$arcilla)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5982 -0.1471 -0.0296  0.1008  0.4786 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.0165     0.3699  10.858 8.64e-09 ***
## ARC$arcilla  -3.6129     1.2228  -2.955  0.00932 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2748 on 16 degrees of freedom
## Multiple R-squared:  0.353,  Adjusted R-squared:  0.3126 
## F-statistic:  8.73 on 1 and 16 DF,  p-value: 0.00932

sshres=shapiro.test(G$residuals)

sshres=shapiro.test(residuales)
ifelse(shres$p.value<0.05, "Rechazo Ho","No rechazo Ho")

## [1] "No rechazo Ho"

Homogeneidad de la varianza.

ARC$Trt <- interaction(ARC$Dosis_insecticida, ARC$Numero_aplicaciones)
bartlett.test(ancova2$residuals~ARC$Trt)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  ancova2$residuals by ARC$Trt
## Bartlett's K-squared = 5.8077, df = 8, p-value = 0.6688

No se cumple el supuesto puesto que el p-valor 0.02329<0.05.

Independencia de los residuales

plot(G$residuals)

Gráfico de interacción

Dosis_insecticida2= as.factor(ARC$Dosis_insecticida)
Aplicaciones2 = as.factor(ARC$Numero_aplicaciones)
Biomasa2 = as.vector(ARC$biomasa)
D_2 <- data.frame(Dosis_insecticida2, Aplicaciones2, Biomasa2)
str(D_2)

## 'data.frame':    18 obs. of  3 variables:
##  $ Dosis_insecticida2: Factor w/ 3 levels "3.25","3.75",..: 2 1 1 1 1 3 2 2 2 2 ...
##  $ Aplicaciones2     : Factor w/ 3 levels "4","5","6": 1 1 3 2 2 2 3 3 2 2 ...
##  $ Biomasa2          : num  2.71 2.77 2.14 2.56 2.71 ...

library(dplyr)
library(ggplot2)
D_2 %>%
  group_by(Dosis_insecticida2, Aplicaciones2) %>%
  summarise(medias_groups = mean(Biomasa2)) -> MEDIAS_MD

## `summarise()` regrouping output by 'Dosis_insecticida2' (override with `.groups` argument)

ggplot(MEDIAS_MD)+
  aes(x=Dosis_insecticida2,y=medias_groups,col = Aplicaciones2) +
  geom_line(aes(group = Aplicaciones2))+
  geom_point()

Finalmente se obtiene un grafico de interacción casi igual al anterior con interraciones de cada una de las aplicaciónes.

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.

library(readxl)
punto5 <- read_excel("C:/Users/Equipo/Downloads/Tabla punto 5 (1).xlsx")
View(punto5) # Se importan los datos reordenados

df = data.frame(punto5) # Se renombra en un data.frame
head(df)

##   finca Variedad Parcela Test Respuesta
## 1     1        1       1    1      9.76
## 2     1        1       2    1     10.65
## 3     1        1       3    1      6.50
## 4     1        1       4    1      8.08
## 5     1        1       5    1      7.84
## 6     1        1       6    1      9.00

library(daewr)

## Registered S3 method overwritten by 'DoE.base':
##   method           from       
##   factorize.factor conf.design

 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.282   0.608  0.438
## Parcela:finca           1    3.1   3.077   0.354  0.554
## Parcela:finca:Variedad  1    5.2   5.179   0.596  0.443
## Residuals              76  660.7   8.693

#Se usa el codigo propuesto
library(lme4)

## Loading required package: Matrix

## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car

## 
## Attaching package: 'lme4'

## The following object is masked from 'package:daewr':
## 
##     cake

modr3 <- lmer( Respuesta ~ 1 + (1|Parcela) + (1|Parcela:finca)+ (1|Parcela:finca:Variedad), data = df)

## boundary (singular) fit: see ?isSingular

summary(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.92753 -0.39932  0.00922  0.43797  1.65397 
## 
## Random effects:
##  Groups                 Name        Variance Std.Dev.
##  Parcela:finca:Variedad (Intercept) 1.2305   1.1093  
##  Parcela:finca          (Intercept) 0.0000   0.0000  
##  Parcela                (Intercept) 7.0127   2.6481  
##  Residual                           0.8795   0.9378  
## Number of obs: 80, groups:  
## Parcela:finca:Variedad, 60; Parcela:finca, 40; Parcela, 20
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   8.2369     0.6188   13.31
## convergence code: 0
## boundary (singular) fit: see ?isSingular

Calculamos la varianza total y se compara con la varianza debida a parcela y variedad

Varparcela= (100*7.0127)/(7.0127+1.2305+0.8795)
Varparcela

## [1] 76.87088

Varvariedad= (100*1.2305)/(1.2305+7.0127+0.8795)
Varvariedad

## [1] 13.48833

Varparcvarie=Varparcela+Varvariedad
Varparcvarie

## [1] 90.35921

76.87 % de la variación total se debe a la variabilidad entre parcelas, y el 13,54% se debe a la variabilidad de las variedades dentro de la finca mientras que dentro de las parcelas (parcelas:finca) la variabilidad es nula.

Se propone analizar unos gráficos y un modelo de correlación para clarificar o ejemplificar la variabilidad debida a parcela y la variable respuesta.

### Gráficos de comportamiento de las medidas

boxplot(df$Respuesta~df$Parcela,  main="Comportamiento de respuesta en función de la parcela",size=0.05, col='aquamarine3')

Tal como lo muestra el analisis de varianzas existe una gran variabilidad entre las parceleas del diseño y debería evaluarse como ocurre con la interacción finca:variedad:parcela.

cor.test(df$Respuesta,df$Parcela,method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  df$Respuesta and df$Parcela
## t = -0.78479, df = 78, p-value = 0.435
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3023487  0.1338074
## sample estimates:
##         cor 
## -0.08851176

plot(df$Respuesta,df$Parcela,pch=16,col="salmon", main = "Correlación de respuesta y parcela")
abline(df$Respuesta,df$Parcela)

Existe una correlación negativa muy baja ( casi nula) entre la respuesta y la parcela. Los datos del gráfico muestran una linea de ajuste que no es muy fiel al posible modelo. La relación lineal entre los datos es dudosa y se aprecia una gran dispersión.

En el enlace https://cran.r-project.org/web/packages/asbio/asbio.pdf se tienen unos datos de potasio de muestras de suelos medidas en 8 diferentes laboratorios.

Compare descriptivamente (medidas, tablas y gráficos) para representar los datos.

library(asbio)

## Loading required package: tcltk

data(K)
klab<-K
klab

##      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

library(psych)

## 
## Attaching package: 'psych'

## The following object is masked from 'package:asbio':
## 
##     skew

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

describe(klab)

##      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.27

describeBy(klab, klab$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.0

Para interpretar se observaram solo los valores de media, mediana y desviación estándar

meanklab=tapply(klab$K, klab$lab, mean)
sdklab=tapply(klab$K, klab$lab, sd)
medianklab=tapply(klab$K, klab$lab, median);
dfklab=data.frame(meanklab,sdklab, medianklab)

library(DT)

datatable(dfklab, class = 'cell-border stripe',filter = 'top', options = list(
  pageLength = 8, autoWidth = TRUE))

Usando un gráfico de violin

library(ggplot2)
ggplot(data = klab , aes(x = lab , y = K)) + 
  geom_jitter(size = 1, color = 'gray2', alpha = 0.5) +
  geom_violin(aes(fill = lab), color = 'black', alpha = .8 ) +
  geom_boxplot(color = 'black', alpha = 0.7) + 
    stat_summary(fun=mean, geom="point", shape=23, size=2, color = "red")+
  xlab('Laboratorio') + 
  ylab('K (Potasio)') +
  ggtitle('Potasio (K) medido en diferentes laboratorios') + 
  theme(legend.position = "none")

## Warning: Computation failed in `stat_summary()`:
## Can't convert a double vector to function

El gráfico de violin muestra la distribución (dispersión) de los datos y si existen o no atípicos, el punto rojo en los boxplot representa las medias y la linea media las medianas. El laboratorio H presenta mediciones atípicas y su concentración es raznoblemnte pareja, si se usa la media el valor más alto cae en el lab J, si se considera la mediana es mejor E, pero J presenta menor dispersión de datos respecto a E por lo que debería ser más confiable. El G presenta menos variabilidad de sus datos. . Hata aquí se pueden consigerar los laboratorios G y J como los mejores. Sin embargo se debe ver más a fondo con otras pruebas.

¿qué prueba me recomendaría para comparar la medida que usted seleccione? Proponga una solución. Sabiendo que son muestras mezcladas de una misma finca, ¿ Se perciben diferencias en las medidas como consecuencia probable de los laboratorios?

Prueba 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).

\[H_0: La~variable~respuesta~es~la~misma~para~todos~los~laboratorios \\ H_a: Al~menos~un~laboratorio~presenta~una~variable~respuesta~distinta\]

kruskalklab<-kruskal.test(klab$K, klab$lab)
kruskalklab

## 
##  Kruskal-Wallis rank sum test
## 
## data:  klab$K and klab$lab
## Kruskal-Wallis chi-squared = 24.482, df = 7, p-value = 0.000937

ifelse(kruskalklab$p.value<0.05, "Rechazo Ho", "No rechazo Ho")

## [1] "Rechazo Ho"

qchisq(0.05, 7, lower.tail = F) # Hallando el valor máximo aceptable para no rechazar la hipotesis nula

## [1] 14.06714

kruskalklab$statistic<qchisq(0.05, 7, lower.tail = F)# Comparando el resultado estadístico de Kruskal con el valor

## Kruskal-Wallis chi-squared 
##                      FALSE

Se rechaza la igualdad en la medición de los laboratorios con un p.valor muy bajo y menor que 0.05 y de los gráficos usados puede deberse principalmente a la medición de H y para evaluar esa diferencia realizamos:

Ahora usando la prueba de Comparación por pares de Nemenyi

PMCMR::posthoc.kruskal.nemenyi.test(klab$K~klab$lab) # se usa la comparación por pares para ver si exusten diferencias entre las mediciones de potasio en los laboratorios

## 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:  klab$K by klab$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

Conclusion

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 laboratios no muestran diferencias significativas estadísticamente hablando sustentados por las ´pruebas de Kruskal y la de Nemenyi, y que es bien representada por el gráfico de violin.

Diseñe un experimento en parcelas divididas en bloques completos (diseño en franjas o strip plot design). Genere los datos usted mismo y esquematice el diseño. Expliqué las razones de colocar el primer factor en la parcela principal y el segundo en la subparcela. Genere unos datos asociados a una covariable y corra el análisis de covarianza respectivo.¿ se justifica el uso de la covariable en el modelo? ¿se justifica el bloque en el modelo? ¿ se tiene interacción de factores? De no presentarse interacción , reduzca el modelo a la presencia de solo términos cuyos p_ valores sean menores al 6%. Escriba el modelo final e interprete el resultado desde un punto de vista agronómico seleccionando el “mejor tratamiento” en la mejor condición de bloqueo y con la presencia de la covariable. No olvide ordenar datos de la covariable. Revise los supuestos necesarios para el análisis estadístico que está proponiendo.

https://psfaculty.plantsciences.ucdavis.edu/agr205/Lectures/2011_Lectures/L12a_SplitPlot.pdf) http://pages.cs.wisc.edu/~songwang/disc11.pdf https://stat.ethz.ch/education/semesters/as2015/anova/08_Split_Plots.pdf