Estimación de la media global \(\mu\)

Si \(\hat{\mu}\) es un estimador para la media global \(\mu\) entonces:

\[\hat{\mu}=\bar{y_{..}}\]

Estimación para el efecto del \(i-ésimo\) tratamiento \(\tau_i\)

Si \(\hat{\tau_i}\) es un estimador para el efecto \(\tau\) del tratamiento \(i\) entonces

\[\hat{\tau_i}=\bar{y_{i.}}-\bar{y_{..}} \\ i= 1,2,...,a\]

Estimación para el efecto del \(j-ésimo\) bloque \(\beta_j\)

Si \(\hat{\beta_j}\) es un estimador para el efecto \(\beta\) del bloque \(j\) entonces:

\[\hat{\beta_j}=\bar{y_{.j}}-\bar{y_{..}}\\ j = 1,2,...,b\]

Estimación de la observación \(y_{ij}\)

Si \(\hat{y_{ij}}\) es un estimador para la observación \(y_{ij}\), entonces

\[\hat{y_{ij}}= \hat{\mu} + \hat{\tau_i} + \hat{\beta_j} \]

Remplazando los estimados de \(\hat{\mu}\), \(\hat{\tau_i}\) y \(\hat{\beta_j}\) obtenemos que:

\[\hat{y_{ij}} = \bar{y_{i.}} + \bar{y_{.j}} + \bar{y_{..}}\]

Verificación gráfica

# ---------- Realizamos el análisis de varianza --------------
library(readxl)
datos <- read_excel("datos.xlsx")
datos$tratamiento <- as.factor(datos$tratamiento)
datos$bloque <- as.factor(datos$bloque)
datos$bloque <- as.factor(datos$bloque)
modelo <- lm(datos$observacion~(datos$tratamiento+datos$bloque))
anova <- aov (modelo)

# -- Realizamos la gráfica cuantil - cuantil de los residuales del ANOVA -- 
library(car)
qqPlot (anova, main= "Gráfico Cuantil - Cuantil", id=FALSE)

Del gráfico Cuantil-Cuantil no se observa alguna evidencia para sospechar que no exista normalidad en los residuales.

Verificación analítica mediante Shapiro-Wilk.

Para la prueba tenemos las siguientes hipótesis:

\[H_0: e_{ij} \sim N(0, \sigma^2) \\ H_1: e_{ij} \nsim N(0, \sigma^2) \\ i,j =1,2,3,4 \]

# ---------- Realizamos el análisis de varianza --------------
library(readxl)
datos <- read_excel("datos.xlsx")
datos$tratamiento <- as.factor(datos$tratamiento)
datos$bloque <- as.factor(datos$bloque)
datos$bloque <- as.factor(datos$bloque)
modelo <- lm(datos$observacion~(datos$tratamiento+datos$bloque))
anova <- aov (modelo)

# -- Realizamos la prueba de Normalidad Shapiro - Wilk -- 
shapiro.test(modelo$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo$residuals
## W = 0.97299, p-value = 0.8844

El estadístico de prueba \(W\) para la prueba de Shapiro - Wilk resultó ser igual a \(W = 0.97299\) se debe comparar con el estadístico teórico \(W_{1-\alpha}\) con nivel de significancia \(\alpha=0.05\) y \(16\) residuales. Este valor de tabla resulta ser \(W_{0.95}=0.981\). Cómo \(W<W{1-\alpha}\) no existe evidencia estadística suficiente para rechazar \(H_o\) por lo que se concluye que los residuales provienen de una distribución normal según el contraste de Shapiro - Wilk

Verificación analítica mediante Kolmogorov - Smirnov.

Para la prueba tenemos las siguientes hipótesis:

\[H_0: e_{ij} \sim N(0, \sigma^2) \\ H_1: e_{ij} \nsim N(0, \sigma^2) \\ i,j =1,2,3,4 \]

# ---------- Realizamos el análisis de varianza --------------
library(readxl)
datos <- read_excel("datos.xlsx")
datos$tratamiento <- as.factor(datos$tratamiento)
datos$bloque <- as.factor(datos$bloque)
datos$bloque <- as.factor(datos$bloque)
modelo <- lm(datos$observacion~(datos$tratamiento+datos$bloque))
anova <- aov (modelo)

# -- Realizamos la prueba de Normalidad Kolmogorov - Smirnov -- 
library(nortest)
lillie.test(modelo$residuals)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  modelo$residuals
## D = 0.11346, p-value = 0.8412

El estadístico de prueba \(D\) para la prueba de Kolmogorov-Smirnov resulta en \(D= 0.11346\) , se compara con el estadístico teórico \(KS\), que resulta, para un nivel de significancia \(\alpha=0.05\) y \(16\) residuales, en:

\[KS = \frac{C_{\alpha}}{k(n)}= \frac{0.895}{4.2025}=0.2129\]

Como \(D<KS\) no se tiene evidencia estadística suficiente para rechazar \(H_0\) por lo que se concluye que los residuales provienen de una distribución normal según el contraste de Kolmogorov - Smirnov.

Verificación gráfica - homocedasticidad.

# ---------- Realizamos el análisis de varianza --------------
library(readxl)
datos <- read_excel("datos.xlsx")
datos$tratamiento <- as.factor(datos$tratamiento)
datos$bloque <- as.factor(datos$bloque)
datos$bloque <- as.factor(datos$bloque)
modelo <- lm(datos$observacion~(datos$tratamiento+datos$bloque))
anova <- aov (modelo)

# ----- Graficamos valores predichos vs Residuales ----- 
predichos <- modelo$fitted.values
residuales <- modelo$residuals
plot(predichos,residuales, main="Predichos VS Residuales", xlab="Valores predichos", ylab="Residuales", ylim = c(-3,3), xlim=c(5,16))

En la verificación gráfica de la homocedasticidad podemos observar cierto patrón tipo “cono”, por lo que resulta importante verificar la homocedasticidad de manera analítica.

Verificación analítica mediante Bartlett

Para la prueba tenemos las siguientes hipótesis:

\[H_0: \sigma_1^2 = \sigma_2^2 = \sigma_3^2 = \sigma_4^2 =\sigma^2 \\ H_1: \sigma_i^2 \neq \sigma_j^2~para~al~menos~un~par~(i,j)~con~i \neq j \\ i,j=1,2,3,4\]

# ---------- Realizamos el análisis de varianza --------------
library(readxl)
datos <- read_excel("datos.xlsx")
datos$tratamiento <- as.factor(datos$tratamiento)
datos$bloque <- as.factor(datos$bloque)
datos$bloque <- as.factor(datos$bloque)
modelo <- lm(datos$observacion~(datos$tratamiento+datos$bloque))
anova <- aov (modelo)

# ----- Prueba de Bartlett ----- 
bartlett.test(anova$residuals~datos$tratamiento)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  anova$residuals by datos$tratamiento
## Bartlett's K-squared = 2.3476, df = 3, p-value = 0.5035

qchisq(0.05,3,lower.tail = FALSE)

## [1] 7.814728

La prueba de Bartlett arrojó un estadístico de prueba \(X_0^2= 2.3476\), el estadístico teórico para nivel de significancia \(\alpha=0.05\) y \(a-1\) grados de libertad de los tratamientos resulta en \(X_{0.05,~3}^2=7.814728\). Como \(X_0^2<X_{\alpha,~a-1}\) entonces no existe evidencia suficiente para rechazar \(H_0\) por lo que las varianzas de los tratamientos son iguales.

Verificación gráfica - independencia.

Para la verificación gráfica de la independencia supondremos que los las observaciones se realizaron en el mismo orden mostrado en la tabla típica de la Situación inicial.

# ---------- Realizamos el análisis de varianza --------------
library(readxl)
datos <- read_excel("datos.xlsx")
datos$tratamiento <- as.factor(datos$tratamiento)
datos$bloque <- as.factor(datos$bloque)
datos$bloque <- as.factor(datos$bloque)
modelo <- lm(datos$observacion~(datos$tratamiento+datos$bloque))
anova <- aov (modelo)

# ----- Grafico residuales como serie de tiempo ----- 
Orden_de_corrida <- c(1:16)
residuales <- anova$residuals
plot(Orden_de_corrida, residuales, xlab="Tiempo", ylab="Residuales", main="Orden de Corrida vs Residuales", xlim = c(0,17), ylim=c(-3,3))

En el gráfico de Orden de corrida VS Residuales no se observa un patrón claro en los datos, por lo que no existe evidencia para dudar sobre la independencia entre los residuales.

Verificación análitica independencia - Prueba Durbin - Watson

Para la realización de la prueba de independencia se tienen las siguientes hipótesis:

\[H_0: \rho = 0 \\ H_1: \rho \neq 0\]

Siendo \(\rho\) el parámetro que representa la correlación entre residuos consecutivos, es decir, \(Correlación(e_t,~e_{t+1})\)

# ---------- Realizamos el análisis de varianza --------------
library(readxl)
datos <- read_excel("datos.xlsx")
datos$tratamiento <- as.factor(datos$tratamiento)
datos$bloque <- as.factor(datos$bloque)
datos$bloque <- as.factor(datos$bloque)
modelo <- lm(datos$observacion~(datos$tratamiento+datos$bloque))
anova <- aov (modelo)

# ----- Prueba de Durbin - Watson ----- 
library(car)
durbinWatsonTest(modelo)

##  lag Autocorrelation D-W Statistic p-value
##    1      -0.4201389      2.711806   0.482
##  Alternative hypothesis: rho != 0

El estadístico de prueba \(d\) para la prueba \(Durbin - Watson\) resultó en \(d=2.711806\), por lo que para un nivel de significancia \(\alpha = 0.05\), \(4~tratamientos\) y \(16~residuales\) y \(d_L=0.734\), \(d_U=1.935\), \(4-d_U=2.065\), \(4-d_L=3.266\). La prueba resulta no ser concluyente.