Parcial Disenio de Experimentos

Valores U = 0, T = 0

Punto 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 803 y que la suma de cuadrados de los tratamientos (SCtrt) es 6000 Con esta información complete la tabla del ANOVA.

\[n=72\\ tratamientos~(trt)=9\\ repeticiones~(rpp)=\frac{n}{trt}=\frac{72}{9}=8\\ Varianza~(S2)=803\\ SC_{trt}=6000\\ Grados~de~libertad~de~los~tratamientos~(gl_{trt})=trt−1=9−1=8\\ Gradosde~libertad~del~error~(gl_{error})=trt*(rpp−1)=9(8−1)=63\\ Grados~de~libertad~totales~(gl_{totales})=n−1=72−1=71\]

Calculando la SC error para completar el ANOVA

\[S^2=\frac{SC_{totales}}{gl}\] Donde: \[SC=Suma~de~cuadrados\] \[S2=Varianza~total\] \[gl=Grados~de~libertad (n−1)\]

Como la suma de cuadrados totales es:

\[SC_{totales}=SC_{trt} + SC_{error}\] Entonces si se reemplaza en la formula de varianza tenemos:

\[S^2=\frac{SC_{trt} + SC_{error}}{gl}\] Donde podemos despejar la suma de cuadrados del error:

\[SC_{error}=(S^2∗gl) −SC_{trt})\]

trt=9
rp=8
n=72
S2=803
SCtrt=6000
gl_trt=8
gl_err=63
gl_tot=71

SCerr=(803*71)-6000;SCerr

## [1] 51013

Es decir que tenemos una sumatoria de cuadrados de error de:

\[SC_{error}=51013\] A partir de las sumatorias de cuadrados ya podemos calcular los cuadrados medios:

CMtrt<- SCtrt/gl_trt;CMtrt

## [1] 750

CMerr<- SCerr/gl_err;CMerr

## [1] 809.7302

\[CM_{trt}=750\] \[CM_{error}=809.73\]

Y posteriormente se calcula el cociente de F:

CF <- CMtrt/CMerr;CF

## [1] 0.9262345

\[Valor~de~F=0.9262345\]

library(readxl)

## Warning: package 'readxl' was built under R version 4.0.3

ANOVA_punto1<- read_excel("D:/Informacion Unal/Semestre 2020-2/Diseño/ANOVA punto1.xlsx")
head(ANOVA_punto1)

## # A tibble: 3 x 5
##   FV         SC    df    CM `Coef F`
##   <chr>   <dbl> <dbl> <dbl>    <dbl>
## 1 Between  6000     8  750     0.901
## 2 Within  52433    63  832.   NA    
## 3 Total   51013    71   NA    NA

Interpretación del valor de F. Con este coeficiente F podemos decir que la variabilidad entre las repeticiones es mayor que la variabilidad entre tratamientos, debido a que su valor es menor a uno.

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 como el cociente F calculado de la tabla para concluir)?

Calculamos el valor de P para la tabla ANOVA:

pf(q = 0.9262,df1 = 8,df2 = 63,lower.tail = F)

## [1] 0.5013187

\[P_{valor} = 0.5013187\]

Hipótesis nula:

\[H_0: \mu_{1}=\mu_{2}=\mu_{3}=\cdots = \mu_{9}\\ H_a: Almenos\ una\ diferente\\Donde~\mu~representa~cada ~tratamiento\]

Realizando el gráfico de la Función densida F de Fisher.

ftab<- 2.8
x <- seq( -4, 4, by = 0.1)
y <- dnorm( x )
plot(function(x) df( x, df1 = 8, df2= 63), 0, 5, ylim = c( 0, 1 ),
      col = "black", type = "l", lwd = 2,
      main = "Función de densidad F de Fisher")
abline(v=ftab,col="yellow")
abline(v=CF, col = "red")
text(2.8,0.4,"Ftab")
text(0.9,0.4,"CF")
text(3.5, 0.6, "Zona de rechazo")
text(1.9, 0.6, "Zona de no rechazo")
text(2.9,0.1,"5%")
text(2.3,0.85,"pvalor = 0.5013")
segments(x0 = 0.9262,y0 = 0.8,x1 = 5, y1 = 0.8, col="green")

Interpretando: Con la gráfica anterior podemos notar que nuestro coeficiente F (CF) se ubica en la región de no rechazo, pues los datos reflejan que estadísticamente son iguales los tratamientos, por eso la hipótesis no se rechaza, con el p valor podríamos llegar a la misma conclusión, aunque estaría tentado a decir con seguridad que los tratamientos son iguales hay que aclarar que esto es un error, debido a que como ya se había mencionado anteriormente, el valor de f es menor a 1 por ende esta igualdad se debe a que hay la variación está dada mayormente por las repeticiones que por los mismos tratamientos.

¿Vale la pena comparar las medias de tratamientos a posteriori del ANOVA (prueba de Tukey)? No, debido a que estadisticamente los tratamientos son iguales, no tendría sentido mirar de forma mas rigurosa la diferencia entre los tratamientos.

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.

Tabla Algodón

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)

library(readxl)
# Se crean factores para realizar el anova
t_algodon <- read_excel("D:/Informacion Unal/Semestre 2020-2/Diseño/Algodon.xlsx")
t_algodon$Granjero = as.factor(t_algodon$Granjero)
t_algodon$Metodo = as.factor(t_algodon$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.design

Utilizamos la función lm para este tipo de diseño.

modelo1 = lm( Perdida ~ Granjero * Metodo, data = t_algodon)
anova(modelo1)

## Warning in anova.lm(modelo1): 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)
## Granjero         5 138.30 27.6608               
## Metodo           4  49.12 12.2799               
## Granjero:Metodo 19  26.23  1.3805               
## Residuals        0   0.00

table (t_algodon$Granjero, t_algodon$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  1

Debido a que sólo hay una respuesta por bloque y una por factor no se puede usar este modelo, pues estipula que es perfecto cuando no lo es.

Para el modelo des utilizaremos la función lm con bloque:

modelo2 <- lm( Perdida ~ Granjero + Metodo, data = t_algodon )
anova( modelo2 )

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

Éste es el método apropiado debido a que se coloca primero el bloque que son los granjeros, ya que nuestra variable de interés es el método.

Si cambiamos el orden entre los bloques y la variable tenemos:

modelo3 <- lm( Perdida ~ Metodo + Granjero, data = t_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

Si se comparan los dos anova anteriores podemos inferir que el orden sí importa, pues aunque ambos muestran valores estadísticos, su interpretación es muy diferente.

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

\[Estimando~el~valor faltante:\\Promedio=\frac{6.75 + 13.05 + 10.26 + 8.01 + 8.42}{5}= 9.298\] Se incluyo el valor x1 en la base de datos, para ello se crea una base de datos con los mismos datos:

t_algodon2 = t_algodon
t_algodon2$Perdida[1] = 9.298

Se calcula el ANOVA de los datos balanceados.

modelo4 = lm( Perdida ~ Granjero + Metodo, data = t_algodon2 )
anova(modelo4)

## Analysis of Variance Table
## 
## Response: Perdida
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Granjero   5 139.364 27.8728 21.0147  2.45e-07 ***
## Metodo     4  51.150 12.7874  9.6411 0.0001637 ***
## Residuals 20  26.527  1.3263                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparando los modelos 2 y 4, podemos notar que hay algunos cambios, especialmente en el valor de F y el P valor, sobre todo el modelo 4 nos arroja un p valor menor, por lo que es preferible

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

set.seed(1911)
#Se fija la semilla para obtener datos 
co10 = runif(50, min = 2, max = 2.2)
co10 = sort.int(co10,partial = 25)
#se ordenan los datos con la funcion sort.int
co5 = runif(50,min =3.0, max = 3.1)
co5 = sort.int(co5,partial = 25)
#se ordenan los datos con la funcion sort.int

window = expand.grid(long = seq(0,100,25), lat = seq(0,200,length.out = 10))
View(window)

codf = data.frame(long = rep(window$long, 2), 
                  lat = rep(window$lat, 2), 
                  prof = rep(c(5, 10), each = 50), 
                  co = c(co5, co10))
View(codf)

library(plotly)

## Warning: package 'plotly' was built under R version 4.0.3

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 4.0.3

## 
## 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 = codf$long, y = codf$lat, z = codf$prof, type = "scatter3d", mode = "markers", color = codf$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.

Para comparar las medias de las capas aplicamos la prueba Paired t-test, con la función t.test, de la siguiente manera:

Hipótesis:

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

Pair_T<-t.test(codf$co~codf$prof,
       alternative='t',
       paired=T)
Pair_T

## 
##  Paired t-test
## 
## data:  codf$co by codf$prof
## t = 189.25, df = 49, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.9348847 0.9549519
## sample estimates:
## mean of the differences 
##               0.9449183

ifelse(Pair_T$p.value<0.05,'medias diferentes','medias iguales')

## [1] "medias diferentes"

De acuerdo con esta prueba podemos notar que las medias de contenido de carbono orgánico entre las capas es diferente.

Punto 4

El siguiente diseño se corresponde con un factorial completo (3^2) en arreglo completamente al azar. Los factores y la respuesta fueron creados con el código:

F1: Dosis de un insecticida que se supone genera efecto de disminución de crecimiento.

F2: Número de aplicaciones durante el desarrollo del cultivo.

D<-expand.grid(F1=c(3.25,3.75,4.25),F2=c(4,5,6)) # Se crea el Diseño 3^2
D<-rbind(D, D)
set.seed(2020)
D<-D[order(sample(1:18)),] # Se aleatoriza la estructura
class(D)

## [1] "data.frame"

D$biomasa=sort.int(rnorm(18,3,0.3),partial = 9)
D

##      F1 F2  biomasa
## 2  3.75  4 2.708826
## 10 3.25  4 2.772692
## 16 3.25  6 2.143359
## 4  3.25  5 2.560519
## 13 3.25  5 2.708666
## 6  4.25  5 2.773705
## 17 3.75  6 2.770350
## 8  3.75  6 2.832470
## 14 3.75  5 2.898280
## 5  3.75  5 3.359619
## 9  4.25  6 3.054099
## 1  3.25  4 3.157896
## 12 4.25  4 3.487669
## 11 3.75  4 3.451547
## 15 4.25  5 3.016111
## 7  3.25  6 3.042156
## 3  4.25  4 3.200552
## 18 4.25  6 2.989329

• Escriba (completamente especificado) el modelo del diseño.

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

• Use los resultados del ANOVA y el gráfico de interacción (ggplot2) para visualizar si existe o no interacción entre los factores.

\[Y_{ijk}=\mu+\alpha_i+\beta_j+(\alpha\beta){ij}+\epsilon_{ijk}\\ i=1,2\\ j=1,2,3\\ k=1,2,\cdots,18\]

D$F1<-as.factor(D$F1); D$F1

##  [1] 3.75 3.25 3.25 3.25 3.25 4.25 3.75 3.75 3.75 3.75 4.25 3.25 4.25 3.75 4.25
## [16] 3.25 4.25 4.25
## Levels: 3.25 3.75 4.25

D$F2<-as.factor(D$F2); D$F2

##  [1] 4 4 6 5 5 5 6 6 5 5 6 4 4 4 5 6 4 6
## Levels: 4 5 6

mod1<- aov(biomasa~F1*F2,data = D);mod1

## Call:
##    aov(formula = biomasa ~ F1 * F2, data = D)
## 
## Terms:
##                        F1        F2     F1:F2 Residuals
## Sum of Squares  0.4160823 0.3425590 0.1634908 0.9459429
## Deg. of Freedom         2         2         4         9
## 
## Residual standard error: 0.3241987
## Estimated effects may be unbalanced

anova(mod1)

## Analysis of Variance Table
## 
## Response: biomasa
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## F1         2 0.41608 0.208041  1.9794 0.1939
## F2         2 0.34256 0.171280  1.6296 0.2489
## F1:F2      4 0.16349 0.040873  0.3889 0.8116
## Residuals  9 0.94594 0.105105

library(ggplot2)
library(gplots)

## Warning: package 'gplots' was built under R version 4.0.3

## 
## Attaching package: 'gplots'

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

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.0.3

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

# Summarizing data
D02 <- D %>% 
  group_by(F1, F2) %>% 
  summarise(biom_mean = mean(biomasa))

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

#y_se = psych::describe(y)$se)

D02 %>% 
  ggplot(aes(F1,biom_mean,color=F2))+
  geom_line(aes(group = F2))

plot(TukeyHSD(mod1))

#El cero está dentro de los tres intervalos de confianza

D$arcilla<-runif(18, min = 0.20, max = 0.40)
table(D$arcilla)

## 
## 0.211400281731039 0.233473422378302 0.245705667557195 0.257527577877045 
##                 1                 1                 1                 1 
## 0.259058687230572 0.277518581692129 0.284879770968109 0.294059629784897 
##                 1                 1                 1                 1 
## 0.306059261178598 0.311170172505081 0.315368745615706 0.322743925265968 
##                 1                 1                 1                 1 
## 0.322866390412673 0.330207164213061  0.33452609446831 0.372663450846449 
##                 1                 1                 1                 1 
## 0.376887316023931 0.398512358311564 
##                 1                 1

aovarc<-aov(biomasa~F1*F2+arcilla,data = D)
summary(aovarc)

##             Df Sum Sq Mean Sq F value Pr(>F)
## F1           2 0.4161 0.20804   1.762  0.232
## F2           2 0.3426 0.17128   1.451  0.290
## arcilla      1 0.0125 0.01253   0.106  0.753
## F1:F2        4 0.1525 0.03813   0.323  0.855
## Residuals    8 0.9444 0.11805

\[Y_{ijk}=\mu+\alpha_i+\beta_j+\ t_{auk}+(\alpha\beta){ij}+\epsilon_{ijk}\\ i=1,2,3\\ j=1,2,3\\ k=1,2,\cdots,18\]

library(multcomp)

## Warning: package 'multcomp' was built under R version 4.0.3

## Loading required package: mvtnorm

## Warning: package 'mvtnorm' was built under R version 4.0.3

## Loading required package: survival

## Loading required package: TH.data

## Warning: package 'TH.data' was built under R version 4.0.3

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:plotly':
## 
##     select

## The following objects are masked from 'package:daewr':
## 
##     cement, chem

## 
## Attaching package: 'TH.data'

## The following object is masked from 'package:MASS':
## 
##     geyser

ancova<- aov(biomasa~F1*F2+arcilla,data = D)
ancova

## Call:
##    aov(formula = biomasa ~ F1 * F2 + arcilla, data = D)
## 
## Terms:
##                        F1        F2   arcilla     F1:F2 Residuals
## Sum of Squares  0.4160823 0.3425590 0.0125326 0.1525359 0.9443652
## Deg. of Freedom         2         2         1         4         8
## 
## Residual standard error: 0.3435777
## Estimated effects may be unbalanced

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.

Diseño Anidado Escalonado

library(readxl)
C_papa <- read_excel("D:/Informacion Unal/Semestre 2020-2/Diseño/Cultivo_de_papas.xlsx")
View(C_papa)
head(C_papa)

## # A tibble: 6 x 5
##   Finca Variedad  Test Parcela Respuesta
##   <dbl>    <dbl> <dbl>   <dbl> <chr>    
## 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

dfpapa = data.frame(C_papa)

library(daewr)
modelo_papa2<-aov(Respuesta ~ Parcela + Parcela:Finca + Parcela:Finca:Variedad, data =
dfpapa)
summary(modelo_papa2)

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

library(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':
## 
##     cake

dfpapa$Respuesta = as.numeric(as.character(dfpapa$Respuesta))

modelo_papalme4 = lmer(Respuesta ~ 1 + (1|Parcela) + (1|Parcela:Finca) +
                + (1|Parcela:Finca:Variedad),data = dfpapa)

## boundary (singular) fit: see ?isSingular

summary(modelo_papalme4)

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## Respuesta ~ 1 + (1 | Parcela) + (1 | Parcela:Finca) + +(1 | Parcela:Finca:Variedad)
##    Data: dfpapa
## 
## 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
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

\[\%_{Varinza~Total}=\frac{(100∗7.0122)}{(7.0122+1.2309+0.8789))}=77\%\] Según estos resultados, el 77% de la variación total se debe exclusivamente a la variabiliada entre parcelas. En cuanto a la interacción entre la parcela y la finca se evidencia que hay una varianza nula.

Diseños de la Librearia Agricolae

library(agricolae)

## Warning: package 'agricolae' was built under R version 4.0.3

library(daewr)
library(lmerTest)

## Warning: package 'lmerTest' was built under R version 4.0.3

## 
## Attaching package: 'lmerTest'

## The following object is masked from 'package:lme4':
## 
##     lmer

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

library(DT)

## Warning: package 'DT' was built under R version 4.0.3

library(ggplot2)

Diseño Grecolatino (Graeco-latin square design)

La función de este diseño es “design.graeco ()”

\[y_{íjkl}=\mu+\theta_i+\tau_j+\omega_k+\psi_l+\epsilon{ijkl}\\ i=1,2,\cdots,\rho\\ j=1,2,\cdots,\rho\\ k=1,2,\cdots,\rho\\ l=1,2,\cdots,\rho\]

Este es un diseño Factorial completo, aleatorizado, sin anidamiento, con bloqueo

Donde: \(y_{íjkl}\) = Observación en fila \(i\), la columna \(l\), para la letra latina \(j\) y la letra griega \(k\) \(\\mu\) = Media global \(\theta_i\) = Efecto de la \(i-ésima\) fila \(\tau_j\) = Efecto del tratamiento \(j\) de la letra latina \(\omega_k\) = Efecto del tratamiento de la letra griega \(k\) \(\psi_l\) = Es el efecto de la colunma \(l\) \(\epsilon_{ijkl}\) = Error aleatorio

temp <- c("T1","T1","T1","T1",
           "T2","T2","T2","T2",
           "T3","T3","T3","T3",
           "T4","T4","T4","T4")
pro <- c("P1","P2","P3","P4",
           "P1","P2","P3","P4",
           "P1","P2","P3","P4",
           "P1","P2","P3","P4")
latin <- c(3,2,4,1,
            2,3,1,4,
            1,4,2,3,
            4,1,3,2)
grec <- c(2,3,4,1,
           1,4,3,2,
           4,1,2,3,
           3,2,1,4)
yi <- c(5,12,13,13,
         6,10,15,11,
         7,5,5,7,
         11,10,8,9)

temp <- factor(temp)
pro <- factor(pro)
latin <- factor(latin)
grec <- factor(grec)

data_1 <- data.frame(temp,pro,latin,grec,yi)
dg <- design.graeco(latin,grec,serie = 0)

## not implemented design  16 x 16 , see help(design.graeco)

graeco_1 <- dg$book
plots <- as.numeric(graeco_1[,1]);plots

## numeric(0)

print(matrix(plots,byrow=TRUE,ncol=4))

##      [,1] [,2] [,3] [,4]

cathe <- lm(yi~pro+temp+latin+grec)
ANOVA <- aov(cathe)
summary(ANOVA)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## pro          3  22.19   7.396   6.017 0.0873 .
## temp         3  57.69  19.229  15.644 0.0245 *
## latin        3  36.69  12.229   9.949 0.0456 *
## grec         3  32.19  10.729   8.729 0.0542 .
## Residuals    3   3.69   1.229                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#pro = proceso #temp = temperatura #latin = letra latina #grec = letra griega

Diseño Strip-plot desing

Strip-Plot desing, con la función “Strip.plot()”

Este diseño es Factorial completo, sin anidamiento, totalmente aleatorizado sin bloqueo

Strip Plot

\[y_{ijk}=\mu+\tau_i+\beta_j+(\tau\beta)_{ij}+\gamma_k+(\tau\gamma)_{ik}+(\beta\gamma)_{jk}+\epsilon_{ijk}\]

library(readxl)
StripPdata <- read_excel("D:/Informacion Unal/Semestre 2020-2/Diseño/StripPlotdata.xlsx")
View(StripPdata)
str(StripPdata)

## tibble [48 x 4] (S3: tbl_df/tbl/data.frame)
##  $ block     : num [1:48] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Varieties : chr [1:48] "V1" "V1" "V1" "V1" ...
##  $ Fertilizer: chr [1:48] "F1" "F2" "F3" "F4" ...
##  $ yield     : num [1:48] 10.2 11.1 6.8 5.3 8 9.7 8.6 3.4 2 10.9 ...

StripPdata$Varieties <- as.factor(StripPdata$Varieties)
StripPdata$Fertilizer <- as.factor(StripPdata$Fertilizer)

str(StripPdata)

## tibble [48 x 4] (S3: tbl_df/tbl/data.frame)
##  $ block     : num [1:48] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Varieties : Factor w/ 3 levels "V1","V2","V3": 1 1 1 1 2 2 2 2 3 3 ...
##  $ Fertilizer: Factor w/ 4 levels "F1","F2","F3",..: 1 2 3 4 1 2 3 4 1 2 ...
##  $ yield     : num [1:48] 10.2 11.1 6.8 5.3 8 9.7 8.6 3.4 2 10.9 ...

attach(StripPdata)

Var <- StripPdata$Varieties
Fert <- StripPdata$Fertilizer
Bloq <- StripPdata$block
Rend <- StripPdata$yield

modSPD = strip.plot(BLOCK = Bloq,
                   COL = Var,
                   ROW = Fert,
                   Y = Rend)

## 
## ANALYSIS STRIP PLOT:  Rend 
## Class level information
## 
## Var  :  V1 V2 V3 
## Fert     :  F1 F2 F3 F4 
## Bloq     :  1 2 3 4 
## 
## Number of observations:  48 
## 
## model Y: Rend ~ Bloq + Var + Ea + Fert + Eb + Fert:Var + Ec 
## 
## Analysis of Variance Table
## 
## Response: Rend
##          Df  Sum Sq Mean Sq F value    Pr(>F)    
## Bloq      3  13.692   4.564  2.4086 0.1007122    
## Var       2 163.007  81.503 57.0397 0.0001248 ***
## Ea        6   8.573   1.429  0.7541 0.6144958    
## Fert      3 152.685  50.895 17.1086 0.0004638 ***
## Eb        9  26.773   2.975  1.5700 0.1983665    
## Fert:Var  6  40.320   6.720  3.5465 0.0170071 *  
## Ec       18  34.107   1.895                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## cv(a) = 16 %, cv(b) = 23 %, cv(c) = 18.4 %, Mean = 7.491667

Punto 8: Resumenes

El uso de Diseños en parcelas divididas:

Es considerado cuando se tienen algunos factores que tienen niveles más difíciles de cambiar que otros, estos se reinician menos seguido que el resto de manera intencional, en un experimento completamente aleatorizado todos los factores se reinician el mismo número de veces. Originados principalmente para configuraciones de agricultura, referidos a parcelas de suelo, factores difíciles se denominan whole-plot, factores fáciles se denominan subplots. Ejemplo, parcelas whole-plot, semillas subplot. #whole-plots = #reinicios; #experimentos = #subplots Son más complicados que un diseño completamente aleatorizado, pero con ventajas estadísticas y de implementación. Esto puede carrear confusiones entre un split-plot inadvertidas al poder ser analizado como un CRD (completely randomized design) que puede llevar a conclusiones erróneas.

Si no se quiere reiniciar todos los factores entre runs, es necesario estructurar mucho mejor el proceso de aleatorización, sin omitirlo; y analizando los datos como un split-plot, esto lleva a conclusiones validas que dan información precisa sobre los efectos. Un análisis correcto demuestra que las observaciones obtenidas sin reinicios en los wp están más relacionadas entre sí que las que involucran reiniciar todos los factores. Para determinar el orden de las runs existen dos aleatorizaciones: el orden para correr los WP y para recolectar información dentro de cada WP. El diseño de parcelas divididas inadvertidas genera un diseño con solo 50% de eficiencia como la mejor opción para agrupar las observaciones dentro de los wp, dejarlo a la suerte es una aproximación con un alto potencial para decepción. Por ende, muchas veces es más eficiente aplicar el split-plot correcto que corriendo el mismo experimento como un CRD.

Los costos son otro factor a considerar, si los factores difíciles de cambiar son caros o requieren de mucho tiempo, escoger un split-plot con menos wp pueden resultar en ahorros importantes; la aproximación a la cuantificación de los costos relativos en un experimento se asume a ser proporcional a la suma del número de wp y el total de observaciones (N) (#WP +rN) r = CSP/CWP C = costo relativo de cambio **dentro del CSP se incluye el costo de la medición y lectura de cada observación. Cuando el factor difícil de cambiar resulta extremadamente más caro que reiniciar el sp (r=0), añadir más observaciones saldría barato y el costo total se vería principalmente afectado por el #reinicios del wp = #WP. r = 0.1-1 costo está dado por ambos factores

Escoger el diseño dependiendo del costo relativo y/o combinar precio y calidad no siempre es apropiado, pero puede conducir a comparaciones más realistas y facilitar la toma de decisiones. Una de las principales ventajas de este tipo de diseños es que los parámetros de ña media del modelo pueden ser estimados por separado de los errores de los WP y los Sp, acción que suele ser difícil de estimar con precisión en otras instancias. Los diseños de parcelas divididas son herramientas importantes y practicas; cuando se escogen estratégicamente, estos diseños pueden incrementar la cantidad de información que se puede extraer de dicho experimento diseñado.

Unidad experimental y unidad de observación:

En la elaboración de diseños experimentales se debe tener mucho cuidado con su estructuración para que se pueda tener una interpretación correcta de los análisis, para ello resulta importante tener en claro lo que es “unidad experimental” y “unidad de observación” La literatura líder en el diseño de experimentos define la unidad experimental, también llamada unidad de replicación, como la entidad más pequeña que se asigna independientemente de todas las demás unidades a un tratamiento particular; la palabra “independiente” es clave para esta definición y una unidad de observación, también conocida como unidad de muestreo, se define como la entidad física en la que se mide un resultado de interés en un experimento. En muchos diseños simples, las unidades experimentales y las unidades de observación son sinónimos; es decir, pueden coincidir con la misma entidad física. Pero en el caso contrario las unidades de observación son dependientes a la unidad experimental, a esto se le llama anidación, cuando esto ocurre las unidades de observación se conocen comúnmente como submuestras, pseudoreplicaciones o réplicas técnicas, para indicar que estas observaciones están correlacionadas y, por lo tanto, no constituyen una verdadera replicación independiente. Se podría decir que las unidades experimentales se definen en términos de asignaciones de tratamiento independientes, mientras que las unidades de observación se definen en términos de mediciones de resultados. Estos son criterios de definición claramente diferentes. Como tal, las observaciones no representan necesariamente replicaciones. Sin embargo, las unidades de observación se encuentran generalmente dentro de unidades experimentales, que a su vez determinan la cantidad de replicación de un experimento dado.

La conclusión es que el modelo estadístico debe describir un proceso plausible que dé lugar a las observaciones mediante (1) capturar las variables independientes importantes que afectan a la variable de resultado, y por (2) especificar cualquier restricción en la aleatorización o cualquier otra estructura de datos que induzca la correlación entre las observaciones. Debe quedar claro que es la variación estimada entre las unidades experimentales independientes la que determina el error experimental adecuado para evaluar los efectos del tratamiento.

La identificación adecuada de la unidad experimental en un estudio determinado debe basarse en relación de la pregunta de investigación específica —¿Cuál es el efecto de un tratamiento en un resultado específico de interés? — en combinación con los muchos matices logísticos del proceso de recopilación de datos, es decir, el diseño experimental y la estructura de datos correspondiente.

Lo más importante es que vale la pena subrayar que la especificación de la unidad experimental en un experimento determinado no es una cuestión de opinión; más bien, está determinado por cómo se creó el experimento, cómo se recopilaron los datos y el alcance de inferencia previsto. Como tal no se puede definir tan fácilmente la asignación de unidad experimental, sino que depende de cada caso que se desea estudiar, haciendo una buena asignación de estas unidades garantiza que tengamos unos resultados estadísticamente significativos y que sean verídicos.

Punto 9: Análisis de un articulo, seleccionar un artículo científico de una revista de agronomía donde se haya utilizado un diseño en parcelas divididas.

EVALUACIÓN AGRONÓMICA Y ECONÓMICA DE ARVEJA ARBUSTIVA (Pisum sativum L.) EN DIFERENTES ÉPOCAS DE SIEMBRA Y SISTEMAS DE TUTORADO

Hacer las críticas constructivas sobre:

La mención de la estructura factorial

No presentan el modelo de organización, pero se infiere que es estructura factorial incompleta

La razón de colocar cada factor en la parcela principal o en la subparcela

Definen muy bien los factores involucrados (época, tutorado, genotipo, temperatura y precipitación), pero en el arreglo de parcelas se definió el tutorado como parcela completa y los genotipos como sub-parcelas, ya que el tutorado es un factor más caro y complicado de cambiar que el resto y a partir de este se pueden definir el resto de factores, mientras que, por otro lado cambiar el genotipo resulta más fácil y sencillo.

La revisión de supuestos para el análisis de varianza

No exponen la revisión de supuestos pero la construcción del análisis de varianza está muy bien construida

La tabla del análisis de varianza

Nombrada Tabla 1. Esta es muy completa, presenta todos los datos de los factores muy bien organizados y especificados, fácil de leer y comprender, expone los datos con diferencias significativas y datos con diferencias ALTAMENTE significativas

El uso de muchos análisis de varianzas en lugar de uno solo multivariante

Se utilizó un único análisis de varianza que encerraba todas las variables a considerar para el primer análisis, esto es preferible ya que al hacer muchos análisis por separado para cada variable esto puede presentar problemas de superposición de los resultados, lo que significa que solo podría analizar los resultados finales de estos más no la interacción que tienen los diferentes factores entre sí.

El método de comparaciones de medias después del Anova

Posteriormente de la realización del anova se realizó una prueba de Tukey (P<0.05) teniendo en cuanta solamente las medias de ciertas variables como lo fueron los Días de cosecha (DC), la altura de la planta (AP) y el número de vainas por planta (NVP) y se analizó la diferencia honesta significativa.

La interacción de factores

Esta se realizó bloqueando el factor Tutorado y al ser un diseño anidado, esto permitió analizar la interacción del bloque con los sub-bloques, estos siendo genotipo, precipitación y temperatura.

Tutorado/ (Genotipo: Precipitación: Temperatura)

La presencia de bloques

El bloque designado fue el Tutorado y a partir de este se desarrollaron los sub-bloques descritos anteriormente i. El balanceo o desbalanceo No hubo desbalanceo ya que todos los datos estaban completos, por lo cual esto no fue un inconveniente al realizar el análisis estadístico

La definición clara de la unidad experimental

Este es uno de los aspectos a resaltar y es que definieron y expusieron claramente cuál era y cómo se realizó la unidad experimental, esta se constituía por cuatro surcos de 3m de largo, a una distancia de 0,5m entre surco y 0,10m entre sitio, colocando una semilla por sitio y dejando calle de 1m, cada cuatro surcos, para facilitar el manejo agronómico, obteniendo una densidad de 160.000 plantas x ha-1.

Software utilizado y librería específica (en caso de ser R)

En el artículo no se especificó el software utilizado o si quiera si se utilizó alguno.

Otros aspectos que considere de interés

Al final se realizó un análisis del porcentaje de vainas sanas, además se concluyó que la implementación de un tutorado se puede traducir en un incremento del 20-22% en el rendimiento comparado con un sistema sin tutorado.

Y finalmente hacen un análisis económico en el que a primera vista el sistema sin tutorado parece la mejor opción, pero a largo plazo la inversión en los sistemas de tutorado se recupera debido al incremento significativo en el rendimiento, a pesar de que el sistema que en un principio generó más ingresos netos fue el sistema sin tutorado.