Experimento de camarones I
Kenneth Roy Cabrera Torres
Jueves, 5 de abril de 2018
- Lectura y adecuación de la base de datos
- Análisis exploratorio
- Modelización
- Modelo final
- Gráficas diagnósticas y verificación de supuestos
- Verificación de supuestos
- Inferencia del modelo
- Inferencia numérica
- Modelo heterocedástico (mínimos cuadrados ponderados)
- Diagnósticos del modelo mínimos cuadrados ponderados
- Inferencia del modelo
Se llevó a cabo un diseño experimental que consistió en tres réplicas para cada una de las 12 combinaciones de tres factores. Los factores que se evaluaron fueron Temperatura (25 ºC y 35 ºC), Densidad de población (80 indv/m\(^3\) y 160 indv/m\(^3\)) y salinidad del agua (10%, 25% y 40%). Se asignaron 36 acuarios al azar en los 36 experimentos que se realizaron en total durante cuatro semanas y se anotaron el aumento de peso en miligramos por camarón cosechado.
Fuente: Diseño de experimentos, Robert O. Kuehl, Thomsom Editores 2001. Tabla 6.13 página 201
#------------------------------------------------------------------------------#
# 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")
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 ...Análisis exploratorio
with(camarones, mean(Aumento.de.peso))## [1] 277.2222with(camarones, sd(Aumento.de.peso))## [1] 125.9086with(camarones, tapply(Aumento.de.peso,
list(Salinidad, Densidad.de.población, Temperatura),
mean))## , , 25 °C
##
## 80 idv/m3 160 idv/m3
## 10% 70.33333 70.66667
## 25% 465.66667 333.00000
## 40% 359.00000 252.33333
##
## , , 35 °C
##
## 80 idv/m3 160 idv/m3
## 10% 408.0000 331.0000
## 25% 274.6667 311.6667
## 40% 243.0000 207.3333with(camarones, tapply(Aumento.de.peso,
list(Salinidad, Densidad.de.población, Temperatura),
sd))## , , 25 °C
##
## 80 idv/m3 160 idv/m3
## 10% 17.15615 16.62328
## 25% 87.64892 108.28204
## 40% 59.85817 11.37248
##
## , , 35 °C
##
## 80 idv/m3 160 idv/m3
## 10% 51.11751 30.11644
## 25% 47.96179 42.66536
## 40% 36.16628 82.62163require(ggplot2)## Loading required package: ggplot2colores <- c("blue","red")
g1 <- ggplot(camarones, aes(Salinidad, Aumento.de.peso, colour = Temperatura))
g1 + geom_point() +
facet_grid(. ~ Densidad.de.población) +
stat_summary(fun.y = mean, geom = "line", aes(group = Temperatura)) +
scale_colour_manual(values = colores) +
ylab("Aumento de peso (mg)") +
ggtitle("Crecimiento de camarones") 
g2 <- ggplot(camarones, aes(Salinidad, Aumento.de.peso, colour = Densidad.de.población))
g2 + geom_point() +
facet_grid(. ~ Temperatura) +
stat_summary(fun.y = mean, geom = "line", aes(group = Densidad.de.población)) +
scale_colour_manual(values = colores) +
ylab("Aumento de peso (mg)") +
ggtitle("Crecimiento de camarones") 
Modelización
El primer modelo propuesto es:
\[ y_{ijkl} = \mu + \alpha_i + \beta_j + \gamma_k + (\alpha \beta)_{ij} + (\alpha \gamma)_{ik} + (\beta \gamma)_{jk} + (\alpha \beta \gamma)_{ijk} + \epsilon_{ijkl} \]
Donde: \[ \epsilon_{ijkl} \sim \mathcal{N}(0, \sigma^2) \quad i.i.d. \]
m1<-aov(Aumento.de.peso ~ Temperatura + Salinidad + Densidad.de.población +
Salinidad:Densidad.de.población +
Temperatura:Salinidad +
Temperatura:Densidad.de.población +
Temperatura:Salinidad:Densidad.de.población,
data=camarones)
summary(m1)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperatura 1 12619 12619 3.837 0.0619 .
## Salinidad 2 98143 49072 14.921 6.15e-05 ***
## Densidad.de.población 1 24754 24754 7.527 0.0113 *
## Salinidad:Densidad.de.población 2 1713 856 0.260 0.7729
## Temperatura:Salinidad 2 308839 154420 46.953 5.06e-09 ***
## Temperatura:Densidad.de.población 1 6669 6669 2.028 0.1673
## Temperatura:Salinidad:Densidad.de.población 2 23187 11593 3.525 0.0455 *
## Residuals 24 78931 3289
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m2<-aov(Aumento.de.peso ~ Temperatura + Salinidad + Densidad.de.población +
Temperatura:Salinidad +
Temperatura:Densidad.de.población +
Temperatura:Salinidad:Densidad.de.población,
data=camarones)
summary(m2)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperatura 1 12619 12619 3.837 0.0619 .
## Salinidad 2 98143 49072 14.921 6.15e-05 ***
## Densidad.de.población 1 24754 24754 7.527 0.0113 *
## Temperatura:Salinidad 2 308839 154420 46.953 5.06e-09 ***
## Temperatura:Densidad.de.población 1 6669 6669 2.028 0.1673
## Temperatura:Salinidad:Densidad.de.población 4 24899 6225 1.893 0.1444
## Residuals 24 78931 3289
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m3<-aov(Aumento.de.peso ~ Temperatura + Salinidad + Densidad.de.población +
Temperatura:Salinidad +
Temperatura:Densidad.de.población,
data=camarones)
summary(m3)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperatura 1 12619 12619 3.403 0.0757 .
## Salinidad 2 98143 49072 13.233 9.00e-05 ***
## Densidad.de.población 1 24754 24754 6.675 0.0153 *
## Temperatura:Salinidad 2 308839 154420 41.643 4.07e-09 ***
## Temperatura:Densidad.de.población 1 6669 6669 1.799 0.1907
## Residuals 28 103830 3708
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m4<-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 ' ' 1Modelo final
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)
g1 <- ggplot(camarones, aes(Salinidad, rstandard(m4)))
g1 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estandarizados")
g1 <- ggplot(camarones, aes(Densidad.de.población, rstandard(m4)))
g1 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estandarizados")
g1 <- ggplot(camarones, aes(Temperatura, rstandard(m4)))
g1 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estandarizados")
g2 <- ggplot(camarones, aes(Salinidad, rstandard(m4)))
g2 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estandarizados") +
facet_grid(. ~ Temperatura)
g2 <- ggplot(camarones, aes(Salinidad, rstandard(m4)))
g2 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estandarizados") +
facet_grid(Densidad.de.población ~ Temperatura)
g1 <- ggplot(camarones, aes(predict(m4), rstandard(m4)))
g1 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estandarizados") +
xlab("Valores estimados bajo el modelo")
op <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
plot(m4)
par(op)Verificación de supuestos
# Normalidad de los residuales
# Generalmente se hace sobre los residuales puros.
shapiro.test(residuals(m4))##
## Shapiro-Wilk normality test
##
## data: residuals(m4)
## W = 0.97834, p-value = 0.6893# Para verificar si existen influencia de datos atípicos sobre la normalidad
shapiro.test(rstandard(m4))##
## Shapiro-Wilk normality test
##
## data: rstandard(m4)
## W = 0.97834, p-value = 0.6893hist(rstandard(m4), freq = FALSE,
main ="Histograma de residuales estandarizados",
las = 1,
xlab = "Residuales estandarizados",
ylab = "Densidad")
lines(density(rstandard(m4)), col = "red", lwd = 2)
# Igualdad de varianzas
bartlett.test(rstandard(m4) ~ Salinidad, data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstandard(m4) by Salinidad
## Bartlett's K-squared = 6.9088, df = 2, p-value = 0.03161bartlett.test(rstandard(m4) ~ Temperatura, data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstandard(m4) by Temperatura
## Bartlett's K-squared = 0.59258, df = 1, p-value = 0.4414bartlett.test(rstandard(m4) ~ Densidad.de.población, data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstandard(m4) by Densidad.de.población
## Bartlett's K-squared = 0.082633, df = 1, p-value = 0.7738bartlett.test(rstandard(m4) ~ paste(Salinidad,Temperatura), data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstandard(m4) by paste(Salinidad, Temperatura)
## Bartlett's K-squared = 7.2421, df = 5, p-value = 0.2032# Para verificar si existen influencia de datos atípicos sobre la normalidad
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) ~ Temperatura, data = camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstudent(m4) by Temperatura
## Bartlett's K-squared = 1.1591, df = 1, p-value = 0.2817bartlett.test(rstudent(m4) ~ Densidad.de.población, data = camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstudent(m4) by Densidad.de.población
## Bartlett's K-squared = 0.31171, df = 1, p-value = 0.5766bartlett.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.09901# Independencia
# Prueba de independencia de Durbin-Watson
require(lmtest)## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericdwtest(m4, order.by = ~ Acuario, data = camarones,
alternative="two.sided")##
## Durbin-Watson test
##
## data: m4
## DW = 1.6151, p-value = 0.2067
## alternative hypothesis: true autocorrelation is not 0Inferencia 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 6# 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)Modelo heterocedástico (mínimos cuadrados ponderados)
Suponiendo que cada combinación tiene su propia varianza:
\[ y_{ijkl} = \mu + \alpha_i + \beta_j + \gamma_k + (\alpha \beta)_{ij} + (\alpha \gamma)_{ik} + (\beta \gamma)_{jk} + (\alpha \beta \gamma)_{ijk} + \epsilon_{ijkl} \]
Donde: \[ \epsilon_{ijkl} \sim \mathcal{N}(0, \sigma_{ijk}^2) \quad i.i.d. \]
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, unioncamarones <- camarones %>%
group_by(Salinidad, Densidad.de.población, Temperatura) %>%
mutate(w = 1 / var(Aumento.de.peso))m1w <- aov(Aumento.de.peso ~ Temperatura + Salinidad + Densidad.de.población +
Salinidad:Densidad.de.población +
Temperatura:Salinidad +
Temperatura:Densidad.de.población +
Temperatura:Salinidad:Densidad.de.población,
data=camarones, weights = w)
summary(m1w)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperatura 1 142.89 142.89 142.888 1.35e-11 ***
## Salinidad 2 276.81 138.40 138.403 6.65e-14 ***
## Densidad.de.población 1 6.37 6.37 6.368 0.01865 *
## Salinidad:Densidad.de.población 2 14.57 7.29 7.287 0.00337 **
## Temperatura:Salinidad 2 157.69 78.84 78.843 2.82e-11 ***
## Temperatura:Densidad.de.población 1 0.26 0.26 0.262 0.61310
## Temperatura:Salinidad:Densidad.de.población 2 9.07 4.54 4.535 0.02134 *
## Residuals 24 24.00 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m2w <- aov(Aumento.de.peso ~ Temperatura + Salinidad + Densidad.de.población +
Salinidad:Densidad.de.población +
Temperatura:Salinidad +
Temperatura:Salinidad:Densidad.de.población,
data=camarones, weights = w)
summary(m2w)## Df Sum Sq Mean Sq F value Pr(>F)
## Temperatura 1 142.89 142.89 142.888 1.35e-11 ***
## Salinidad 2 276.81 138.40 138.403 6.65e-14 ***
## Densidad.de.población 1 6.37 6.37 6.368 0.01865 *
## Salinidad:Densidad.de.población 2 14.57 7.29 7.287 0.00337 **
## Temperatura:Salinidad 2 157.69 78.84 78.843 2.82e-11 ***
## Temperatura:Salinidad:Densidad.de.población 3 9.33 3.11 3.111 0.04518 *
## Residuals 24 24.00 1.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Diagnósticos del modelo mínimos cuadrados ponderados
require(ggplot2)
g1 <- ggplot(camarones, aes(Salinidad, rstandard(m2w)))
g1 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estandarizados")
g1 <- ggplot(camarones, aes(Densidad.de.población, rstandard(m2w)))
g1 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estandarizados")
g1 <- ggplot(camarones, aes(Temperatura, rstandard(m2w)))
g1 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estandarizados")
g2 <- ggplot(camarones, aes(Salinidad, rstandard(m2w)))
g2 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estandarizados") +
facet_grid(. ~ Temperatura)
g2 <- ggplot(camarones, aes(Salinidad, rstandard(m2w)))
g2 + geom_point(col = "blue") +
geom_abline(intercept = c(-2, 2), slope = c(0, 0),
col = "blue", linetype = "dashed") +
ylab("Residuales estandarizados") +
facet_grid(Densidad.de.población ~ Temperatura)
g1 <- ggplot(camarones, aes(predict(m2w), rstandard(m2w)))
g1 + geom_point(col = "blue") +
ylab("Residuales estandarizados") +
xlab("Valores estimados bajo el modelo")
require(ggplot2)
g1 <- ggplot(camarones, aes(Salinidad, residuals(m2w)))
g1 + geom_point(col = "blue") +
ylab("Residuales")
g1 <- ggplot(camarones, aes(Densidad.de.población, residuals(m2w)))
g1 + geom_point(col = "blue") +
ylab("Residuales")
g1 <- ggplot(camarones, aes(Temperatura, residuals(m2w)))
g1 + geom_point(col = "blue") +
ylab("Residuales")
g2 <- ggplot(camarones, aes(Salinidad, residuals(m2w)))
g2 + geom_point(col = "blue") +
ylab("Residuales") +
facet_grid(. ~ Temperatura)
g2 <- ggplot(camarones, aes(Salinidad, residuals(m2w)))
g2 + geom_point(col = "blue") +
ylab("Residuales") +
facet_grid(Densidad.de.población ~ Temperatura)
g1 <- ggplot(camarones, aes(predict(m2w), residuals(m2w)))
g1 + geom_point(col = "blue") +
ylab("Residuales") +
xlab("Valores estimados bajo el modelo")
op <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
plot(m2w)
par(op)# Normalidad de los residuales
# Generalmente se hace sobre los residuales puros.
shapiro.test(residuals(m2w))##
## Shapiro-Wilk normality test
##
## data: residuals(m2w)
## W = 0.96264, p-value = 0.259# Sobre los resuduales estandarizados (no se aplica porque está suponiendo igual varianza)
shapiro.test(rstandard(m2w))##
## Shapiro-Wilk normality test
##
## data: rstandard(m2w)
## W = 0.88163, p-value = 0.001114hist(residuals(m2w), freq = FALSE,
main ="Histograma de residuales",
las = 1,
xlab = "Residuales estandarizados",
ylab = "Densidad")
lines(density(residuals(m2w)), col = "red", lwd = 2)
hist(rstandard(m2w), freq = FALSE,
main ="Histograma de residuales estandarizados",
las = 1,
xlab = "Residuales estandarizados",
ylab = "Densidad")
lines(density(rstandard(m2w)), col = "red", lwd = 2)
# Igualdad de varianzas
bartlett.test(residuals(m2w) ~ Salinidad, data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: residuals(m2w) by Salinidad
## Bartlett's K-squared = 7.3109, df = 2, p-value = 0.02585bartlett.test(residuals(m2w) ~ Temperatura, data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: residuals(m2w) by Temperatura
## Bartlett's K-squared = 0.67803, df = 1, p-value = 0.4103bartlett.test(residuals(m2w) ~ Densidad.de.población, data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: residuals(m2w) by Densidad.de.población
## Bartlett's K-squared = 0.16234, df = 1, p-value = 0.687bartlett.test(residuals(m2w) ~ paste(Salinidad,Temperatura), data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: residuals(m2w) by paste(Salinidad, Temperatura)
## Bartlett's K-squared = 13.059, df = 5, p-value = 0.02283# Igualdad de varianzas
bartlett.test(rstandard(m2w) ~ Salinidad, data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstandard(m2w) by Salinidad
## Bartlett's K-squared = -4.2684e-15, df = 2, p-value = 1bartlett.test(rstandard(m2w) ~ Temperatura, data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstandard(m2w) by Temperatura
## Bartlett's K-squared = -2.157e-16, df = 1, p-value = 1bartlett.test(rstandard(m2w) ~ Densidad.de.población, data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstandard(m2w) by Densidad.de.población
## Bartlett's K-squared = -3.4512e-15, df = 1, p-value = 1bartlett.test(rstandard(m2w) ~ paste(Salinidad,Temperatura), data=camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstandard(m2w) by paste(Salinidad, Temperatura)
## Bartlett's K-squared = -2.4722e-15, df = 5, p-value = 1#
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodeleveneTest(rstandard(m2w) ~ Salinidad, data=camarones)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.0037 0.9963
## 33leveneTest(rstandard(m2w) ~ Temperatura, data=camarones)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0122 0.9128
## 34leveneTest(rstandard(m2w) ~ Densidad.de.población, data=camarones)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0021 0.9641
## 34leveneTest(rstandard(m2w) ~ paste(Salinidad,Temperatura), data=camarones)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.041 0.999
## 30# Para verificar si existen influencia de datos atípicos sobre la normalidad
shapiro.test(rstudent(m2w))##
## Shapiro-Wilk normality test
##
## data: rstudent(m2w)
## W = 0.88601, p-value = 0.001446# Igualdad de varianzas
bartlett.test(rstudent(m2w) ~ Salinidad, data = camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstudent(m2w) by Salinidad
## Bartlett's K-squared = 8.6156e-08, df = 2, p-value = 1bartlett.test(rstudent(m2w) ~ Temperatura, data = camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstudent(m2w) by Temperatura
## Bartlett's K-squared = 5.7649e-11, df = 1, p-value = 1bartlett.test(rstudent(m2w) ~ Densidad.de.población, data = camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstudent(m2w) by Densidad.de.población
## Bartlett's K-squared = 1.1e-07, df = 1, p-value = 0.9997bartlett.test(rstudent(m2w) ~ paste(Salinidad, Temperatura), data = camarones)##
## Bartlett test of homogeneity of variances
##
## data: rstudent(m2w) by paste(Salinidad, Temperatura)
## Bartlett's K-squared = 7.3906e-07, df = 5, p-value = 1library(car)
leveneTest(rstudent(m2w) ~ Salinidad, data = camarones)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.0034 0.9966
## 33leveneTest(rstudent(m2w) ~ Temperatura, data = camarones)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0113 0.9159
## 34leveneTest(rstudent(m2w) ~ Densidad.de.población, data = camarones)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.002 0.9647
## 34leveneTest(rstudent(m2w) ~ paste(Salinidad, Temperatura), data = camarones)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.0377 0.9991
## 30# Comparaciones pareadas mediante Tukey
tk <- TukeyHSD(m2w)
# Gráfica de comparaciones pareada de Tukey en la interacción
op <- par(no.readonly = TRUE)
par(mar=op$mar + c(0,14,0,0))
plot(tk,las=1)
par(op)tk <- TukeyHSD(m2w, "Temperatura:Salinidad:Densidad.de.población")
op <- par(no.readonly = TRUE)
par(mar=op$mar + c(0,15,0,0))
plot(tk, las = 1)
par(op)Inferencia del modelo
# Adición de la predicción del modelo a la base de datos
camarones <- data.frame(camarones[, 1:5], predict(m2w, 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)")
library(dplyr)
camarones$residuales <- residuals(m2w)
camRes <- camarones %>%
group_by(Salinidad, Densidad.de.población, Temperatura) %>%
summarise(sdRes = sd(residuales))
camRes## # A tibble: 12 x 4
## # Groups: Salinidad, Densidad.de.población [?]
## Salinidad Densidad.de.población Temperatura sdRes
## <fct> <fct> <fct> <dbl>
## 1 10% 80 idv/m3 25 °C 17.2
## 2 10% 80 idv/m3 35 °C 51.1
## 3 10% 160 idv/m3 25 °C 16.6
## 4 10% 160 idv/m3 35 °C 30.1
## 5 25% 80 idv/m3 25 °C 87.6
## 6 25% 80 idv/m3 35 °C 48.0
## 7 25% 160 idv/m3 25 °C 108.
## 8 25% 160 idv/m3 35 °C 42.7
## 9 40% 80 idv/m3 25 °C 59.9
## 10 40% 80 idv/m3 35 °C 36.2
## 11 40% 160 idv/m3 25 °C 11.4
## 12 40% 160 idv/m3 35 °C 82.6camRes %>% arrange(sdRes)## # A tibble: 12 x 4
## # Groups: Salinidad, Densidad.de.población [6]
## Salinidad Densidad.de.población Temperatura sdRes
## <fct> <fct> <fct> <dbl>
## 1 40% 160 idv/m3 25 °C 11.4
## 2 10% 160 idv/m3 25 °C 16.6
## 3 10% 80 idv/m3 25 °C 17.2
## 4 10% 160 idv/m3 35 °C 30.1
## 5 40% 80 idv/m3 35 °C 36.2
## 6 25% 160 idv/m3 35 °C 42.7
## 7 25% 80 idv/m3 35 °C 48.0
## 8 10% 80 idv/m3 35 °C 51.1
## 9 40% 80 idv/m3 25 °C 59.9
## 10 40% 160 idv/m3 35 °C 82.6
## 11 25% 80 idv/m3 25 °C 87.6
## 12 25% 160 idv/m3 25 °C 108.