Experimento de camarones II
Kenneth Roy Cabrera Torres
Jueves, 5 de abril de 2018
library(multcomp)## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserlibrary(multcompView)
library(sandwich)#------------------------------------------------------------------------------#
# Funciones auxiliares para hallar el error estándar de los resdiuales
# de un modelo "aov"
options(width=120)
errsdred<-function(modelo){
aovm<-anova(modelo)
return(sqrt(aovm$"Mean Sq"[length(aovm$"Mean Sq")]))
}Lectura y adecuación de la base de datos
camarones<-read.csv2("camarones.csv", fileEncoding = "latin1")
str(camarones)## 'data.frame': 36 obs. of 5 variables:
## $ Acuario : int 16 6 25 8 22 18 11 7 31 28 ...
## $ Aumento.de.peso : int 86 52 73 544 371 482 390 290 397 53 ...
## $ Salinidad : int 10 10 10 25 25 25 40 40 40 10 ...
## $ Densidad.de.población: int 80 80 80 80 80 80 80 80 80 160 ...
## $ Temperatura : int 25 25 25 25 25 25 25 25 25 25 ...# Convertir a factores las variables consideradas factores
camarones$Salinidad<-factor(camarones$Salinidad,
labels = c("10%","25%","40%"))
camarones$Densidad.de.población<-factor(camarones$Densidad.de.población,
levels=c("80","160"),
labels = c("80 idv/m3","160 idv/m3"))
camarones$Temperatura<-factor(camarones$Temperatura,
labels = c("25 °C","35 °C"))
str(camarones)## 'data.frame': 36 obs. of 5 variables:
## $ Acuario : int 16 6 25 8 22 18 11 7 31 28 ...
## $ Aumento.de.peso : int 86 52 73 544 371 482 390 290 397 53 ...
## $ Salinidad : Factor w/ 3 levels "10%","25%","40%": 1 1 1 2 2 2 3 3 3 1 ...
## $ Densidad.de.población: Factor w/ 2 levels "80 idv/m3","160 idv/m3": 1 1 1 1 1 1 1 1 1 2 ...
## $ Temperatura : Factor w/ 2 levels "25 °C","35 °C": 1 1 1 1 1 1 1 1 1 1 ...Modelización
El modelo final es:
\[ y_{ijkl} = \mu + \alpha_i + \beta_j + \gamma_k + (\alpha \beta)_{ij} + \epsilon_{ijkl} \]
Donde: \[ \epsilon_{ijkl} \sim \mathcal{N}(0, \sigma_{ij}^2) i.i.d. \]
Tomemos:
\[ \begin{eqnarray*} \alpha_i & = & \text{Temperatura} \\ \beta_j & = & \text{Salinidad} \\ \gamma_k & = & \text{Densidad de población} \\ \end{eqnarray*} \]
m4 <- aov(Aumento.de.peso~Temperatura+Salinidad+Densidad.de.población+
Temperatura:Salinidad,
data=camarones)
summary(m4)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperatura 1 12619 12619 3.312 0.0791 .
## Salinidad 2 98143 49072 12.879 9.94e-05 ***
## Densidad.de.población 1 24754 24754 6.496 0.0164 *
## Temperatura:Salinidad 2 308839 154420 40.527 4.00e-09 ***
## Residuals 29 110500 3810
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Gráficas diagnósticas y verificación de supuestos
require(ggplot2)## Loading required package: ggplot2g1 <- ggplot(camarones, aes(Salinidad,rstudent(m4)))
g1 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estudentizados")
g1 <- ggplot(camarones, aes(Salinidad,rstudent(m4)))
g1 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
facet_grid(. ~ Temperatura)
ylab("Residuales estudentizados")## $y
## [1] "Residuales estudentizados"
##
## attr(,"class")
## [1] "labels"Verificación de supuestos
# Normalidad de los residuales
# Se recomienda hacer sobre los residuales estudentizados.
shapiro.test(rstudent(m4))##
## Shapiro-Wilk normality test
##
## data: rstudent(m4)
## W = 0.96047, p-value = 0.2222# Igualdad de varianzas
bartlett.test(rstudent(m4)~Salinidad,data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstudent(m4) by Salinidad
## Bartlett's K-squared = 8.7048, df = 2, p-value = 0.01288bartlett.test(rstudent(m4)~paste(Salinidad,Temperatura),data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstudent(m4) by paste(Salinidad, Temperatura)
## Bartlett's K-squared = 9.2633, df = 5, p-value = 0.09901Inferencia del modelo
# Adición de la predicción del modelo a la base de datos
camarones <- cbind(camarones, predict(m4, interval = "confidence"))
# Gráficas teniendo en cuenta el modelo final
require(ggplot2)
colores <- c("blue", "red")
g1 <- ggplot(camarones, aes(Salinidad, Aumento.de.peso, colour = Temperatura))
g1 + geom_point() +
facet_grid(. ~ Densidad.de.población) +
geom_line(aes(Salinidad, fit, group = Temperatura)) +
geom_errorbar(aes(ymax = upr, ymin = lwr), width = 0.25, alpha = 0.2) +
stat_summary(aes(Salinidad, fit, group = Temperatura), fun.y = mean,
geom = "point", shape = 17, size = 3) +
scale_colour_manual(values = colores) +
ylab("Aumento de peso (mg)")
g1 <- ggplot(camarones, aes(Salinidad, Aumento.de.peso,
colour = Densidad.de.población))
g1 + geom_point() +
facet_grid(. ~ Temperatura) +
geom_line(aes(Salinidad, fit, group = Densidad.de.población)) +
geom_errorbar(aes(ymax = upr, ymin = lwr), width = 0.25, alpha = 0.2) +
stat_summary(aes(Salinidad, fit, group = Densidad.de.población),
fun.y = mean,
geom = "point", shape = 17, size = 3) +
scale_colour_manual(values = colores) +
ylab("Aumento de peso (mg)")
Inferencia numérica
model.tables(m4, se = TRUE)## Tables of effects
##
## Temperatura
## Temperatura
## 25 °C 35 °C
## -18.722 18.722
##
## Salinidad
## Salinidad
## 10% 25% 40%
## -57.22 69.03 -11.81
##
## Densidad.de.población
## Densidad.de.población
## 80 idv/m3 160 idv/m3
## 26.222 -26.222
##
## Temperatura:Salinidad
## Salinidad
## Temperatura 10% 25% 40%
## 25 °C -130.78 71.81 58.97
## 35 °C 130.78 -71.81 -58.97
##
## Standard errors of effects
## Temperatura Salinidad Densidad.de.población Temperatura:Salinidad
## 14.55 17.82 14.55 25.20
## replic. 18 12 18 6model.tables(m4,type="means", se = TRUE)## Tables of means
## Grand mean
##
## 277.2222
##
## Temperatura
## Temperatura
## 25 °C 35 °C
## 258.50 295.94
##
## Salinidad
## Salinidad
## 10% 25% 40%
## 220.0 346.2 265.4
##
## Densidad.de.población
## Densidad.de.población
## 80 idv/m3 160 idv/m3
## 303.44 251.00
##
## Temperatura:Salinidad
## Salinidad
## Temperatura 10% 25% 40%
## 25 °C 70.5 399.3 305.7
## 35 °C 369.5 293.2 225.2
##
## Standard errors for differences of means
## Temperatura Salinidad Densidad.de.población Temperatura:Salinidad
## 20.58 25.20 20.58 35.64
## replic. 18 12 18 6with(camarones, tapply(fit, list(Salinidad, Temperatura,
Densidad.de.población),
mean))## , , 80 idv/m3
##
## 25 °C 35 °C
## 10% 96.72222 395.7222
## 25% 425.55556 319.3889
## 40% 331.88889 251.3889
##
## , , 160 idv/m3
##
## 25 °C 35 °C
## 10% 44.27778 343.2778
## 25% 373.11111 266.9444
## 40% 279.44444 198.9444with(camarones, tapply(lwr, list(Salinidad, Temperatura,
Densidad.de.población),
mean))## , , 80 idv/m3
##
## 25 °C 35 °C
## 10% 41.05224 340.0522
## 25% 369.88557 263.7189
## 40% 276.21890 195.7189
##
## , , 160 idv/m3
##
## 25 °C 35 °C
## 10% -11.39221 287.6078
## 25% 317.44112 211.2745
## 40% 223.77446 143.2745with(camarones, tapply(upr, list(Salinidad, Temperatura,
Densidad.de.población),
mean))## , , 80 idv/m3
##
## 25 °C 35 °C
## 10% 152.3922 451.3922
## 25% 481.2255 375.0589
## 40% 387.5589 307.0589
##
## , , 160 idv/m3
##
## 25 °C 35 °C
## 10% 99.94776 398.9478
## 25% 428.78110 322.6144
## 40% 335.11443 254.6144# Error de los residuales (Estimación de la varianza del error experimental)
errsdred(m4)## [1] 61.72787# Comparaciones pareadas mediante Tukey
tk<-TukeyHSD(m4)
# Gráfica de comparaciones pareada de Tukey en la interacción
op <- par(no.readonly = TRUE)
par(mar=op$mar+c(0,6,0,0))
plot(tk,las=1)
par(op)Comparaciones pareadas bajo heterocedasticidad.
require(multcomp)
require(multcompView)
require(sandwich)
compMultHetero <- glht(m4, linfct = mcp(Temperatura = "Tukey"),
vcov = sandwich)## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found -- default contrast might be inappropriatesummary(compMultHetero)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = Aumento.de.peso ~ Temperatura + Salinidad + Densidad.de.población +
## Temperatura:Salinidad, data = camarones)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 35 °C - 25 °C == 0 299.00 19.19 15.58 1.33e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)op <- par(no.readonly = TRUE)
par(mar=op$mar+c(0,6,0,0))
plot(compMultHetero)
par(op)require(multcomp)
require(multcompView)
require(sandwich)
compMultHetero <- glht(m4, linfct = mcp(Densidad.de.población = "Tukey"),
vcov = sandwich)
summary(compMultHetero)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = Aumento.de.peso ~ Temperatura + Salinidad + Densidad.de.población +
## Temperatura:Salinidad, data = camarones)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 160 idv/m3 - 80 idv/m3 == 0 -52.44 18.47 -2.84 0.00817 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)op <- par(no.readonly = TRUE)
par(mar=op$mar+c(0,6,0,0))
plot(compMultHetero)
par(op)require(multcomp)
require(multcompView)
require(sandwich)
compMultHetero <- glht(m4, linfct = mcp(Salinidad = "Tukey"),
vcov = sandwich)## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found -- default contrast might be inappropriatesummary(compMultHetero)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = Aumento.de.peso ~ Temperatura + Salinidad + Densidad.de.población +
## Temperatura:Salinidad, data = camarones)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 25% - 10% == 0 328.83 38.65 8.507 <0.001 ***
## 40% - 10% == 0 235.17 21.83 10.773 <0.001 ***
## 40% - 25% == 0 -93.67 40.93 -2.289 0.0694 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)op <- par(no.readonly = TRUE)
par(mar=op$mar+c(0,6,0,0))
plot(compMultHetero)
par(op)require(multcomp)
require(multcompView)
require(sandwich)
# Nueva variable que representa la interacción
camarones$Temp_Sali <- with(camarones, interaction(Temperatura, Salinidad))
m5 <- aov(Aumento.de.peso ~ -1 +
Temp_Sali, data = camarones)
summary(m5)## Df Sum Sq Mean Sq F value Pr(>F)
## Temp_Sali 6 3186279 531046 117.8 <2e-16 ***
## Residuals 30 135253 4508
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1compMultHetero <- glht(m5, linfct = mcp(Temp_Sali = "Tukey"),
vcov = sandwich)
summary(compMultHetero)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = Aumento.de.peso ~ -1 + Temp_Sali, data = camarones)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 35 °C.10% - 25 °C.10% == 0 299.00 21.78 13.729 <0.001 ***
## 25 °C.25% - 25 °C.10% == 0 328.83 42.93 7.659 <0.001 ***
## 35 °C.25% - 25 °C.10% == 0 222.67 17.82 12.493 <0.001 ***
## 25 °C.40% - 25 °C.10% == 0 235.17 26.68 8.813 <0.001 ***
## 35 °C.40% - 25 °C.10% == 0 154.67 23.17 6.677 <0.001 ***
## 25 °C.25% - 35 °C.10% == 0 29.83 47.48 0.628 0.9856
## 35 °C.25% - 35 °C.10% == 0 -76.33 26.99 -2.828 0.0707 .
## 25 °C.40% - 35 °C.10% == 0 -63.83 33.51 -1.905 0.3878
## 35 °C.40% - 35 °C.10% == 0 -144.33 30.78 -4.689 <0.001 ***
## 35 °C.25% - 25 °C.25% == 0 -106.17 45.80 -2.318 0.1971
## 25 °C.40% - 25 °C.25% == 0 -93.67 49.92 -1.876 0.4043
## 35 °C.40% - 25 °C.25% == 0 -174.17 48.13 -3.619 0.0110 *
## 25 °C.40% - 35 °C.25% == 0 12.50 31.09 0.402 0.9982
## 35 °C.40% - 35 °C.25% == 0 -68.00 28.12 -2.418 0.1640
## 35 °C.40% - 25 °C.40% == 0 -80.50 34.43 -2.338 0.1901
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)op <- par(no.readonly = TRUE)
par(mar=op$mar+c(0,6,0,0))
plot(compMultHetero)
par(op)