# DCA con diferente numero de repeticiones
# Ejemplo:4Se quiere estudiar la influencia de cuatro tipos de riego en la producci昼㸳n de ma攼㹤z. Para ello,
#se seleccionan parcelas del mismo tama昼㸱o y con la misma calidad de terreno, y se prueban los cuatro tipos de riego,
#en 6, 8, 9 y 7 parcelas respectivamente. Los datos de producci昼㸳n que se obtienen son los siguientes:
#Tabla 3.8: Producci昼㸳n de ma攼㹤z por parcela del mismo tama昼㸱o y calidad.
#Tipo de riego // Producci昼㸳n
#A 47, 42, 43, 46, 44, 42
#B 51, 58, 62, 49, 53, 51, 50, 49
#C 37, 39, 41, 38, 39, 37, 42, 36, 40
#D 42, 43, 42, 45, 47, 50, 48
#trat=triego
#tiempo=produccion
produccion<-c(47, 42, 43, 46, 44, 42,51, 58, 62, 49, 53, 51, 50, 49,
37, 39, 41, 38, 39, 37, 42, 36, 40,42, 43, 42, 45, 47, 50, 48)
triego<-c(rep(1,6),rep(2,8),rep(3,9),rep(4,7))
triego
## [1] 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4
# creando un data.frame
datos<-data.frame(triego=factor(triego,labels = c("RIEGO A","RIEGO B","RIEGO C","RIEGO D")),produccion)
datos
## triego produccion
## 1 RIEGO A 47
## 2 RIEGO A 42
## 3 RIEGO A 43
## 4 RIEGO A 46
## 5 RIEGO A 44
## 6 RIEGO A 42
## 7 RIEGO B 51
## 8 RIEGO B 58
## 9 RIEGO B 62
## 10 RIEGO B 49
## 11 RIEGO B 53
## 12 RIEGO B 51
## 13 RIEGO B 50
## 14 RIEGO B 49
## 15 RIEGO C 37
## 16 RIEGO C 39
## 17 RIEGO C 41
## 18 RIEGO C 38
## 19 RIEGO C 39
## 20 RIEGO C 37
## 21 RIEGO C 42
## 22 RIEGO C 36
## 23 RIEGO C 40
## 24 RIEGO D 42
## 25 RIEGO D 43
## 26 RIEGO D 42
## 27 RIEGO D 45
## 28 RIEGO D 47
## 29 RIEGO D 50
## 30 RIEGO D 48
# Diagrama de cajas
plot(datos,col=terrain.colors(4))
# cuadro ANVA
resp<-aov(produccion~triego,data=datos)
anova(resp)
## Analysis of Variance Table
##
## Response: produccion
## Df Sum Sq Mean Sq F value Pr(>F)
## triego 3 850.84 283.614 27.529 3.137e-08 ***
## Residuals 26 267.86 10.302
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# residuales e IC
resp<-lm(produccion~triego-1,data=datos)
summary(resp)
##
## Call:
## lm(formula = produccion ~ triego - 1, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8750 -2.0000 -0.5317 1.9286 9.1250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## triegoRIEGO A 44.000 1.310 33.58 <2e-16 ***
## triegoRIEGO B 52.875 1.135 46.59 <2e-16 ***
## triegoRIEGO C 38.778 1.070 36.24 <2e-16 ***
## triegoRIEGO D 45.286 1.213 37.33 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.21 on 26 degrees of freedom
## Multiple R-squared: 0.9957, Adjusted R-squared: 0.995
## F-statistic: 1501 on 4 and 26 DF, p-value: < 2.2e-16
# intervalos de confianza
confint(resp,level=0.99)
## 0.5 % 99.5 %
## triegoRIEGO A 40.35888 47.64112
## triegoRIEGO B 49.72170 56.02830
## triegoRIEGO C 35.80482 41.75074
## triegoRIEGO D 41.91469 48.65674
# histograma de frecuencias
h<-hist(resp$residuals, col=terrain.colors(10))
xfit <- seq(min(resp$residuals), max(resp$residuals),length = 20)
yfit <- dnorm(xfit,mean=mean(resp$residuals),sd=sd(resp$residuals))
yfit <- yfit*diff(h$mids[1:2])*length(resp$residuals)
lines(xfit, yfit, col = "blue", lwd = 2)
# graficos de probabilidad normal
qqnorm(resp$residuals)
qqline(resp$residuals,col="blue", lwd=2)
# test de normalidad
shapiro.test(resp$residuals)
##
## Shapiro-Wilk normality test
##
## data: resp$residuals
## W = 0.91939, p-value = 0.02586
# grafico de residuales
plot(resp$fitted.values,resp$residuals,main="Estimado vs
Residual")
abline(h=0,lty=2)
#
plot(resp$residuals~triego,main="Tratamiento vs Residual")
#test de homogeneidad de varianzas
bartlett.test(produccion~triego)
##
## Bartlett test of homogeneity of variances
##
## data: produccion by triego
## Bartlett's K-squared = 6.3963, df = 3, p-value = 0.09384
# independencia
plot(resp$residuals,xlab="Orden de corrida")
abline(h=0,lty=2)
# test de Durwin Watson
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(lmtest)
dwtest(resp)
##
## Durbin-Watson test
##
## data: resp
## DW = 1.6827, p-value = 0.06935
## alternative hypothesis: true autocorrelation is greater than 0
# comparaciones multiples
library(agricolae)
# LSDL
LSD.test(resp,"trat",alpha=0.05,console=TRUE)
## Name: trat
## triego
# duncan test
duncan.test(resp,"trat",alpha=0.05,console=TRUE)
## Name: trat
## triego
# tukey test
cmT <- HSD.test(resp, "trat", console=TRUE)
## Name: trat
## triego