setwd("~/Google Drive/Agrosavia/Colaboraciones/Lucero/data")
regaliasdata = read.table("muestra4.csv", header=T, sep=",")
attach(regaliasdata)
library(nlme) #activates the nlme library
library(ggplot2)
## pH loess-time-course figure
pht<- ggplot(regaliasdata, aes(x = time.num.x)) +
facet_grid(gen.x~trat.x) +
geom_point(aes(y = ph.testa, colour = "ph.testa")) +
geom_smooth(aes(y = ph.testa, colour = "ph.testa"), method = "loess", se=T) +
scale_y_continuous(name = "pH") + # Etiqueta de la variable continua
scale_x_continuous(name = "day", breaks=seq(0,7,1)) + # Etiqueta de los grupos
ggtitle("pH time course") + # TÃtulo del plot
theme(axis.line = element_line(colour = "black", # Personalización del tema
size = 0.25))
pht + geom_smooth(aes(y = ph.grano, colour = "ph.grano"), method = "loess", se=T) +
geom_point(aes(y=ph.grano, colour = "ph.grano"))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Acidity loess-time-course figure
act1<- ggplot(regaliasdata, aes(x = time.num.x)) +
facet_grid(gen.x~trat.x) +
geom_point(aes(y=acid.testa, colour = "acid.testa")) +
geom_smooth(aes(y = acid.testa, colour = "acid.testa"), method = "loess", se=T) +
scale_y_continuous(name = "Total acidity") + # Etiqueta de la variable continua
scale_x_continuous(name = "day", breaks=seq(0,7,1)) + # Etiqueta de los grupos
ggtitle("Acidity time course") + # TÃtulo del plot
theme(axis.line = element_line(colour = "black", # Personalización del tema
size = 0.25))
act1 + geom_smooth(aes(y = acid.grano, colour = "acid.grano"), method = "loess", se=T) +
geom_point(aes(y=acid.grano, colour = "acid.grano"))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## Temperature differences time course
temp1<- ggplot(regaliasdata, aes(x = time.num.x)) +
facet_grid(gen.x~trat.x) +
geom_point(aes(y=temp, colour = "temp")) +
geom_smooth(aes(y = temp, colour = "temp"), method = "loess", se=T) +
scale_y_continuous(name = "Temperature °C") + # Etiqueta de la variable continua
scale_x_continuous(name = "day", breaks=seq(0,7,1)) + # Etiqueta de los grupos
ggtitle("Temperature time course") + # TÃtulo del plot
theme(axis.line = element_line(colour = "black", # Personalización del tema
size = 0.25))
temp1 + geom_line(aes(y = tt, colour = "tt")) +
geom_point(aes(y=tt, colour = "tt"))
## `geom_smooth()` using formula 'y ~ x'
## Repeated measures analysis
aov.temp <- aov(temp~factor(trat.x)*factor(gen.x)*factor(time.num.x) + Error(factor(muestra.num.x)))
aov.ph.grano <- aov(ph.grano~factor(trat.x)*factor(gen.x)*factor(time.num.x) + Error(factor(muestra.num.x)))
aov.ph.testa <- aov(ph.testa~factor(trat.x)*factor(gen.x)*factor(time.num.x) + Error(factor(muestra.num.x)))
summary(aov.temp)
##
## Error: factor(muestra.num.x)
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(trat.x) 1 9.41 9.41 16.49 0.000323 ***
## factor(gen.x) 2 120.67 60.34 105.73 2.60e-14 ***
## factor(trat.x):factor(gen.x) 2 22.97 11.48 20.12 2.87e-06 ***
## Residuals 30 17.12 0.57
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value
## factor(time.num.x) 7 13153 1879.0 3328.34
## factor(trat.x):factor(time.num.x) 7 195 27.9 49.41
## factor(gen.x):factor(time.num.x) 14 181 12.9 22.93
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 142 10.1 17.97
## Residuals 210 119 0.6
## Pr(>F)
## factor(time.num.x) <2e-16 ***
## factor(trat.x):factor(time.num.x) <2e-16 ***
## factor(gen.x):factor(time.num.x) <2e-16 ***
## factor(trat.x):factor(gen.x):factor(time.num.x) <2e-16 ***
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov.ph.grano)
##
## Error: factor(muestra.num.x)
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(trat.x) 1 24.464 24.464 440.487 < 2e-16 ***
## factor(gen.x) 2 0.525 0.263 4.727 0.0164 *
## factor(trat.x):factor(gen.x) 2 10.555 5.277 95.024 1.04e-13 ***
## Residuals 30 1.666 0.056
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value
## factor(time.num.x) 7 38.65 5.521 111.09
## factor(trat.x):factor(time.num.x) 7 21.31 3.044 61.25
## factor(gen.x):factor(time.num.x) 14 12.78 0.913 18.36
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 13.55 0.968 19.48
## Residuals 210 10.44 0.050
## Pr(>F)
## factor(time.num.x) <2e-16 ***
## factor(trat.x):factor(time.num.x) <2e-16 ***
## factor(gen.x):factor(time.num.x) <2e-16 ***
## factor(trat.x):factor(gen.x):factor(time.num.x) <2e-16 ***
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov.ph.testa)
##
## Error: factor(muestra.num.x)
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(trat.x) 1 39.68 39.68 476.8 < 2e-16 ***
## factor(gen.x) 2 23.40 11.70 140.6 5.79e-16 ***
## factor(trat.x):factor(gen.x) 2 21.99 11.00 132.1 1.34e-15 ***
## Residuals 30 2.50 0.08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value
## factor(time.num.x) 7 334.6 47.80 1120.78
## factor(trat.x):factor(time.num.x) 7 15.6 2.23 52.23
## factor(gen.x):factor(time.num.x) 14 19.7 1.40 32.91
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 20.5 1.47 34.42
## Residuals 210 9.0 0.04
## Pr(>F)
## factor(time.num.x) <2e-16 ***
## factor(trat.x):factor(time.num.x) <2e-16 ***
## factor(gen.x):factor(time.num.x) <2e-16 ***
## factor(trat.x):factor(gen.x):factor(time.num.x) <2e-16 ***
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Temperature Generalized least squares fit by REML and correlations with different structures
interaction.plot (time.num.x, factor(gen.x), temp, lty=c(1:3),lwd=2,ylab="mean of temp", xlab="time", trace.label="Genotype")
interaction.plot (time.num.x, factor(trat.x), temp, lty=c(1:3),lwd=2,ylab="mean of temp", xlab="time", trace.label="treatment")
nestinginfo <- groupedData(temp ~ gen.x | muestra.num.x, data= regaliasdata)
fit.compsym <- gls(temp~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corCompSymm(, form= ~ 1 | muestra.num.x))
fit.nostruct <- gls(temp~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corSymm(, form= ~ 1 | muestra.num.x), weights = varIdent(form = ~ 1 | time.num.x))
fit.ar1 <- gls(temp~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x))
fit.ar1het <- gls(temp~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x), weights=varIdent(form = ~ 1 | time.num.x))
anova(fit.compsym, fit.nostruct, fit.ar1, fit.ar1het) #compares the models
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.compsym 1 50 730.2069 904.2388 -315.1035
## fit.nostruct 2 84 682.5646 974.9383 -257.2823 1 vs 2 115.64232 <.0001
## fit.ar1 3 50 730.0494 904.0814 -315.0247 2 vs 3 115.48487 <.0001
## fit.ar1het 4 57 688.9371 887.3335 -287.4685 3 vs 4 55.11236 <.0001
## Qualitative assessment of model source significance
anova(fit.nostruct)
## Denom. DF: 240
## numDF F-value p-value
## (Intercept) 1 2197906.0 <.0001
## factor(trat.x) 1 254.7 <.0001
## factor(gen.x) 2 356.6 <.0001
## factor(time.num.x) 7 16033.8 <.0001
## factor(trat.x):factor(gen.x) 2 43.2 <.0001
## factor(trat.x):factor(time.num.x) 7 147.3 <.0001
## factor(gen.x):factor(time.num.x) 14 25.1 <.0001
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 44.4 <.0001
anova(fit.ar1het)
## Denom. DF: 240
## numDF F-value p-value
## (Intercept) 1 1082592.9 <.0001
## factor(trat.x) 1 2.0 0.1632
## factor(gen.x) 2 174.9 <.0001
## factor(time.num.x) 7 8402.1 <.0001
## factor(trat.x):factor(gen.x) 2 17.1 <.0001
## factor(trat.x):factor(time.num.x) 7 86.0 <.0001
## factor(gen.x):factor(time.num.x) 14 21.0 <.0001
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 34.1 <.0001
## Final model diagnostic plot
plot(fit.ar1het)
shapiro.test(resid(fit.ar1het))
##
## Shapiro-Wilk normality test
##
## data: resid(fit.ar1het)
## W = 0.87094, p-value = 8.008e-15
qqnorm(fit.ar1het, abline = c(0,1))
plot(fit.ar1het, resid(., type = "p") ~ fitted(.) | muestra.num.x, abline = 0)
plot(fit.ar1het, muestra.num.x ~ resid(.))
plot(fit.ar1het, acid.grano ~ fitted(.) | muestra.num.x, abline = c(0,1))
## pH testa Generalized least squares fit by REML and correlations with different structures
interaction.plot (time.num.x, factor(gen.x), ph.testa, lty=c(1:3),lwd=2,ylab="mean of ph.testa", xlab="time", trace.label="Genotype")
interaction.plot (time.num.x, factor(trat.x), ph.testa, lty=c(1:3),lwd=2,ylab="mean of ph.testa", xlab="time", trace.label="Treatment")
nestinginfo <- groupedData(ph.testa ~ gen.x | muestra.num.x, data= regaliasdata)
fit.compsym <- gls(ph.testa~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corCompSymm(, form= ~ 1 | muestra.num.x))
fit.nostruct <- gls(ph.testa~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corSymm(, form= ~ 1 | muestra.num.x), weights = varIdent(form = ~ 1 | time.num.x))
fit.ar1 <- gls(ph.testa~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x))
fit.ar1het <- gls(ph.testa~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x), weights=varIdent(form = ~ 1 | time.num.x))
anova(fit.compsym, fit.nostruct, fit.ar1, fit.ar1het) #compares the models
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.compsym 1 50 130.0091 304.0411 -15.004556
## fit.nostruct 2 84 112.9344 405.3081 27.532786 1 vs 2 85.07468 <.0001
## fit.ar1 3 50 131.9932 306.0251 -15.996585 2 vs 3 87.05874 <.0001
## fit.ar1het 4 57 110.3315 308.7279 1.834265 3 vs 4 35.66170 <.0001
## Qualitative assessment of model source significance
anova(fit.nostruct)
## Denom. DF: 240
## numDF F-value p-value
## (Intercept) 1 199533.74 <.0001
## factor(trat.x) 1 745.27 <.0001
## factor(gen.x) 2 606.39 <.0001
## factor(time.num.x) 7 585.73 <.0001
## factor(trat.x):factor(gen.x) 2 243.59 <.0001
## factor(trat.x):factor(time.num.x) 7 55.08 <.0001
## factor(gen.x):factor(time.num.x) 14 44.55 <.0001
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 41.06 <.0001
anova(fit.ar1het)
## Denom. DF: 240
## numDF F-value p-value
## (Intercept) 1 156444.83 <.0001
## factor(trat.x) 1 518.92 <.0001
## factor(gen.x) 2 420.51 <.0001
## factor(time.num.x) 7 1035.99 <.0001
## factor(trat.x):factor(gen.x) 2 161.61 <.0001
## factor(trat.x):factor(time.num.x) 7 51.36 <.0001
## factor(gen.x):factor(time.num.x) 14 34.84 <.0001
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 34.81 <.0001
## Final model diagnostic plots
plot(fit.ar1het)
shapiro.test(resid(fit.ar1het))
##
## Shapiro-Wilk normality test
##
## data: resid(fit.ar1het)
## W = 0.92201, p-value = 4.04e-11
qqnorm(fit.ar1het, abline = c(0,1))
plot(fit.ar1het, resid(., type = "p") ~ fitted(.) | muestra.num.x, abline = 0)
plot(fit.ar1het, muestra.num.x ~ resid(.))
plot(fit.ar1het, acid.grano ~ fitted(.) | muestra.num.x, abline = c(0,1))
## pH grano Generalized least squares fit by REML and correlations with different structures
interaction.plot (time.num.x, factor(gen.x), ph.grano, lty=c(1:3),lwd=2,ylab="mean of ph.grano", xlab="time", trace.label="Genotype")
interaction.plot (time.num.x, factor(trat.x), ph.grano, lty=c(1:3),lwd=2,ylab="mean of ph.grano", xlab="time", trace.label="Treatment")
nestinginfo <- groupedData(ph.grano ~ gen.x | muestra.num.x, data= regaliasdata)
fit.compsym <- gls(ph.grano~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corCompSymm(, form= ~ 1 | muestra.num.x))
fit.nostruct <- gls(ph.grano~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corSymm(, form= ~ 1 | muestra.num.x), weights = varIdent(form = ~ 1 | time.num.x))
fit.ar1 <- gls(ph.grano~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x))
fit.ar1het <- gls(ph.grano~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x), weights=varIdent(form = ~ 1 | time.num.x))
anova(fit.compsym, fit.nostruct, fit.ar1, fit.ar1het) #compares the models
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.compsym 1 50 150.00569 324.0376 -25.00284
## fit.nostruct 2 84 79.36411 371.7378 44.31795 1 vs 2 138.64158 <.0001
## fit.ar1 3 50 147.22633 321.2583 -23.61316 2 vs 3 135.86222 <.0001
## fit.ar1het 4 57 107.30664 305.7031 3.34668 3 vs 4 53.91968 <.0001
## Qualitative assessment of model source significance
anova(fit.nostruct)
## Denom. DF: 240
## numDF F-value p-value
## (Intercept) 1 458680.6 <.0001
## factor(trat.x) 1 1151.1 <.0001
## factor(gen.x) 2 4.4 0.0127
## factor(time.num.x) 7 179.9 <.0001
## factor(trat.x):factor(gen.x) 2 244.9 <.0001
## factor(trat.x):factor(time.num.x) 7 99.6 <.0001
## factor(gen.x):factor(time.num.x) 14 14.1 <.0001
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 30.1 <.0001
anova(fit.ar1het)
## Denom. DF: 240
## numDF F-value p-value
## (Intercept) 1 254840.62 <.0001
## factor(trat.x) 1 917.11 <.0001
## factor(gen.x) 2 10.72 <.0001
## factor(time.num.x) 7 178.93 <.0001
## factor(trat.x):factor(gen.x) 2 198.85 <.0001
## factor(trat.x):factor(time.num.x) 7 96.50 <.0001
## factor(gen.x):factor(time.num.x) 14 18.54 <.0001
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 28.79 <.0001
## Final model diagnostic plots
plot(fit.ar1het)
shapiro.test(resid(fit.ar1het))
##
## Shapiro-Wilk normality test
##
## data: resid(fit.ar1het)
## W = 0.87996, p-value = 2.97e-14
qqnorm(fit.ar1het, abline = c(0,1))
plot(fit.ar1het, resid(., type = "p") ~ fitted(.) | muestra.num.x, abline = 0)
plot(fit.ar1het, muestra.num.x ~ resid(.))
plot(fit.ar1het, acid.grano ~ fitted(.) | muestra.num.x, abline = c(0,1))
## Testa acidity Generalized least squares fit by REML and correlations with different structures
interaction.plot (time.num.x, factor(gen.x), acid.testa, lty=c(1:3),lwd=2,ylab="mean of acid.testa", xlab="time", trace.label="Genotype")
interaction.plot (time.num.x, factor(trat.x), acid.testa, lty=c(1:3),lwd=2,ylab="mean of acid.testa", xlab="time", trace.label="Treatment")
nestinginfo <- groupedData(acid.testa ~ gen.x | muestra.num.x, data= regaliasdata)
fit.compsym <- gls(acid.testa~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corCompSymm(, form= ~ 1 | muestra.num.x))
fit.nostruct <- gls(acid.testa~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corSymm(, form= ~ 1 | muestra.num.x), weights = varIdent(form = ~ 1 | time.num.x))
fit.ar1 <- gls(acid.testa~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x))
fit.ar1het <- gls(acid.testa~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x), weights=varIdent(form = ~ 1 | time.num.x))
anova(fit.compsym, fit.nostruct, fit.ar1, fit.ar1het) #compares the models
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.compsym 1 50 -368.1875 -194.1555 234.0937
## fit.nostruct 2 84 -432.4076 -140.0339 300.2038 1 vs 2 132.22010 <.0001
## fit.ar1 3 50 -357.7315 -183.6995 228.8657 2 vs 3 142.67613 <.0001
## fit.ar1het 4 57 -404.7125 -206.3161 259.3563 3 vs 4 60.98109 <.0001
## Qualitative assessment of model source significance
anova(fit.nostruct)
## Denom. DF: 240
## numDF F-value p-value
## (Intercept) 1 14458.726 <.0001
## factor(trat.x) 1 277.713 <.0001
## factor(gen.x) 2 76.355 <.0001
## factor(time.num.x) 7 110.309 <.0001
## factor(trat.x):factor(gen.x) 2 119.085 <.0001
## factor(trat.x):factor(time.num.x) 7 95.639 <.0001
## factor(gen.x):factor(time.num.x) 14 43.479 <.0001
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 43.485 <.0001
anova(fit.ar1het)
## Denom. DF: 240
## numDF F-value p-value
## (Intercept) 1 11403.132 <.0001
## factor(trat.x) 1 153.767 <.0001
## factor(gen.x) 2 81.948 <.0001
## factor(time.num.x) 7 118.615 <.0001
## factor(trat.x):factor(gen.x) 2 70.978 <.0001
## factor(trat.x):factor(time.num.x) 7 119.930 <.0001
## factor(gen.x):factor(time.num.x) 14 34.945 <.0001
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 45.038 <.0001
## Final model diagnostic plots
plot(fit.ar1het)
shapiro.test(resid(fit.ar1het))
##
## Shapiro-Wilk normality test
##
## data: resid(fit.ar1het)
## W = 0.80233, p-value < 2.2e-16
qqnorm(fit.ar1het, abline = c(0,1))
plot(fit.ar1het, resid(., type = "p") ~ fitted(.) | muestra.num.x, abline = 0)
plot(fit.ar1het, muestra.num.x ~ resid(.))
plot(fit.ar1het, acid.grano ~ fitted(.) | muestra.num.x, abline = c(0,1))
## Bean acidity Generalized least squares fit by REML and correlations with different structures
interaction.plot (time.num.x, factor(gen.x), acid.grano, lty=c(1:3),lwd=2,ylab="mean of acid.grano", xlab="time", trace.label="Genotype")
interaction.plot (time.num.x, factor(trat.x), acid.grano, lty=c(1:3),lwd=2,ylab="mean of acid.grano", xlab="time", trace.label="Treatment")
nestinginfo <- groupedData(acid.grano ~ gen.x | muestra.num.x, data= regaliasdata)
fit.compsym <- gls(acid.grano~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corCompSymm(, form= ~ 1 | muestra.num.x))
fit.nostruct <- gls(acid.grano~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corSymm(, form= ~ 1 | muestra.num.x), weights = varIdent(form = ~ 1 | time.num.x))
fit.ar1 <- gls(acid.grano~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x))
fit.ar1het <- gls(acid.grano~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x), weights=varIdent(form = ~ 1 | time.num.x))
anova(fit.compsym, fit.nostruct, fit.ar1, fit.ar1het) #compares the models
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.compsym 1 50 -456.2439 -282.2120 278.1220
## fit.nostruct 2 84 -552.1086 -259.7349 360.0543 1 vs 2 163.86467 <.0001
## fit.ar1 3 50 -452.0177 -277.9857 276.0088 2 vs 3 168.09094 <.0001
## fit.ar1het 4 57 -493.8161 -295.4197 303.9080 3 vs 4 55.79841 <.0001
## Qualitative assessment of model source significance
anova(fit.nostruct)
## Denom. DF: 240
## numDF F-value p-value
## (Intercept) 1 11447.892 <.0001
## factor(trat.x) 1 783.274 <.0001
## factor(gen.x) 2 197.399 <.0001
## factor(time.num.x) 7 455.050 <.0001
## factor(trat.x):factor(gen.x) 2 169.147 <.0001
## factor(trat.x):factor(time.num.x) 7 119.756 <.0001
## factor(gen.x):factor(time.num.x) 14 23.122 <.0001
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 36.011 <.0001
anova(fit.ar1het)
## Denom. DF: 240
## numDF F-value p-value
## (Intercept) 1 4103.024 <.0001
## factor(trat.x) 1 187.478 <.0001
## factor(gen.x) 2 88.618 <.0001
## factor(time.num.x) 7 202.264 <.0001
## factor(trat.x):factor(gen.x) 2 36.058 <.0001
## factor(trat.x):factor(time.num.x) 7 97.298 <.0001
## factor(gen.x):factor(time.num.x) 14 19.549 <.0001
## factor(trat.x):factor(gen.x):factor(time.num.x) 14 30.408 <.0001
## Final model diagnostic plots
plot(fit.ar1het)
shapiro.test(resid(fit.ar1het))
##
## Shapiro-Wilk normality test
##
## data: resid(fit.ar1het)
## W = 0.815, p-value < 2.2e-16
qqnorm(fit.ar1het, abline = c(0,1))
plot(fit.ar1het, resid(., type = "p") ~ fitted(.) | muestra.num.x, abline = 0)
plot(fit.ar1het, muestra.num.x ~ resid(.))
plot(fit.ar1het, acid.grano ~ fitted(.) | muestra.num.x, abline = c(0,1))
## Fermentation profile analysis
regaliasdata4567 <- subset(regaliasdata, regaliasdata$time.num.x == c("4", "5", "6", "7"))
detach(regaliasdata)
attach(regaliasdata4567)
## Fermentation profile Generalized least squares fit by REML and correlations with different structures
interaction.plot (time.num.x, factor(gen.x), ferm, lty=c(1:3),lwd=2,ylab="mean of ferm", xlab="time", trace.label="Genotype")
interaction.plot (time.num.x, factor(trat.x), ferm, lty=c(1:3),lwd=2,ylab="mean of ferm", xlab="time", trace.label="Treatment")
nestinginfo <- groupedData(ferm ~ gen.x | muestra.num.x, data= regaliasdata4567)
fit.nostruct <- gls(ferm~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corSymm(, form= ~ 1 | muestra.num.x), weights = varIdent(form = ~ 1 | time.num.x))
fit.ar1 <- gls(ferm~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x))
fit.ar1het <- gls(ferm~(factor(trat.x)*factor(gen.x))*factor(time.num.x), data=nestinginfo, corr=corAR1(, form= ~ 1 | muestra.num.x), weights=varIdent(form = ~ 1 | time.num.x))
anova(fit.nostruct, fit.ar1, fit.ar1het) #compares the models
## Model df AIC BIC logLik Test L.Ratio p-value
## fit.nostruct 1 28 139.6833 153.2607 -41.84165
## fit.ar1 2 26 138.1802 150.7878 -43.09009 1 vs 2 2.496895 0.2869
## fit.ar1het 3 29 141.6833 155.7456 -41.84165 2 vs 3 2.496895 0.4759
## Qualitative assessment of model source significance
anova(fit.nostruct)
## Denom. DF: 12
## numDF F-value p-value
## (Intercept) 1 5199.403 <.0001
## factor(trat.x) 1 39.603 <.0001
## factor(gen.x) 2 6.992 0.0097
## factor(time.num.x) 3 34.128 <.0001
## factor(trat.x):factor(gen.x) 2 1.607 0.2407
## factor(trat.x):factor(time.num.x) 3 1.134 0.3746
## factor(gen.x):factor(time.num.x) 6 3.358 0.0352
## factor(trat.x):factor(gen.x):factor(time.num.x) 6 0.799 0.5888
anova(fit.ar1het)
## Denom. DF: 12
## numDF F-value p-value
## (Intercept) 1 5199.403 <.0001
## factor(trat.x) 1 39.603 <.0001
## factor(gen.x) 2 6.992 0.0097
## factor(time.num.x) 3 34.128 <.0001
## factor(trat.x):factor(gen.x) 2 1.607 0.2407
## factor(trat.x):factor(time.num.x) 3 1.134 0.3746
## factor(gen.x):factor(time.num.x) 6 3.358 0.0352
## factor(trat.x):factor(gen.x):factor(time.num.x) 6 0.799 0.5888
## Final model diagnostic plots
plot(fit.ar1het)
shapiro.test(resid(fit.ar1het))
##
## Shapiro-Wilk normality test
##
## data: resid(fit.ar1het)
## W = 0.90819, p-value = 0.00576
qqnorm(fit.ar1het, abline = c(0,1))
plot(fit.ar1het, resid(., type = "p") ~ fitted(.) | muestra.num.x, abline = 0)
plot(fit.ar1het, muestra.num.x ~ resid(.))
plot(fit.ar1het, ferm ~ fitted(.) | muestra.num.x, abline = c(0,1))
