RMD Anovas variables fisiologicas
Mariam Vásquez, Pedro Lizarazo
2023-11-06
library(pacman)
p_load(readxl, dplyr, tidyverse, plyr, magrittr, agricolae, tibble,
ggplot2, tidyr, rsample, forestmangr, broom, mice, car, openxlsx, xlsx2dfs)E N S A Y O 2
Cargar base de datos
BD2 <- read_excel("BD2.complete.fisiologicos.xlsx")
BD2.aov <- BD2 %>%
dplyr::select(., -B, -PL, -FC, -R, -CCI, -FC, -CCI.I, -FC.I) %>%
mutate(., TTO=as.factor(TTO), N=as.factor(N), P=as.factor(P), SDT=as.factor(SDT) )
str(BD2.aov)## tibble [805 × 7] (S3: tbl_df/tbl/data.frame)
## $ TTO.SDT: chr [1:805] "0.0N/0.0P_7" "0.0N/0.0P_7" "0.0N/0.0P_7" "0.0N/0.0P_7" ...
## $ TTO : Factor w/ 12 levels "0.0N/0.0P","0.0N/0.25P",..: 1 1 1 1 1 1 1 1 4 4 ...
## $ N : Factor w/ 4 levels "0.0N","1.1N",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ P : Factor w/ 3 levels "0.0P","0.25P",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ SDT : Factor w/ 3 levels "16","25","7": 3 3 3 3 3 3 3 3 3 3 ...
## $ CCI.R : num [1:805] 13.3 13.3 15.4 13.5 11.9 12.5 11.1 10.5 33.5 17.1 ...
## $ FC.R : num [1:805] 0.815 0.833 0.804 0.851 0.729 0.729 0.729 0.729 0.72 0.72 ...Para CCI
Exploratorio CCI
ggplot(BD2, aes(x=factor(SDT, levels = c("7", "16", "25")), y=CCI.R))+
geom_boxplot(alpha=0.2)+
geom_jitter(aes(color=B), width = 0.1, alpha=0.5)+
facet_grid(~TTO) + ylim(c(0, 55))
25 SDT
# Subconjunto base datos
BD2_25 <- as.data.frame(subset(BD2, subset = SDT==25,
select = c(TTO, N, P, CCI.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_CCI2_25 <- data.frame(lam25 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam25 <- seq(from = -1, to = 20, by = 0.5)
for (i in lam25) {
N <- bartlett.test((CCI.R)^i ~ N, data = BD2_25)
P <- bartlett.test((CCI.R)^i ~ P, data = BD2_25)
# Añadir los valores al dataframe
df_CCI2_25 <- rbind(df_CCI2_25, data.frame(lam25 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_CCI2_25$N.p.value~df_CCI2_25$lam25, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_CCI2_25$P.p.value~df_CCI2_25$lam25, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(CCI.R ~ N * P, data = BD2_25, lambda = seq(-1, 20, 1/10))
#lambda.CCI2_25 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.CCI2_25 <- 1
# Creación del modelo ANOVA de dos factores tipo III
mod_CCI.25 <- aov( (CCI.R)^lambda.CCI2_25 ~ N * P, data = BD2_25)
AOV_CCI.25 <- Anova(mod_CCI.25, type = "III"); AOV_CCI.25## Anova Table (Type III tests)
##
## Response: (CCI.R)^lambda.CCI2_25
## Sum Sq Df F value Pr(>F)
## (Intercept) 7228.9 1 141.8776 < 2.2e-16 ***
## N 1877.3 3 12.2812 1.802e-07 ***
## P 244.4 2 2.3983 0.093191 .
## N:P 1069.6 6 3.4988 0.002502 **
## Residuals 11464.2 225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_CCI.25 <- mod_CCI.25$residuals
res.est_CCI.25 <-res_CCI.25/sd(res_CCI.25)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD2_25$CCI.R), main="Variable original" )
hist( (BD2_25$CCI.R)^lambda.CCI2_25, main="Variable transformada" )
boxplot( (BD2_25$CCI.R), horizontal = T, main="Variable original" )
boxplot( (BD2_25$CCI.R)^lambda.CCI2_25, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_CCI.25, which = 2)
library(nortest)
ad.test(res.est_CCI.25)##
## Anderson-Darling normality test
##
## data: res.est_CCI.25
## A = 0.46839, p-value = 0.2472cvm.test(res.est_CCI.25)##
## Cramer-von Mises normality test
##
## data: res.est_CCI.25
## W = 0.070868, p-value = 0.2713pearson.test(res.est_CCI.25)##
## Pearson chi-square normality test
##
## data: res.est_CCI.25
## P = 25.405, p-value = 0.04476ks.test(res.est_CCI.25, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_CCI.25
## D = 0.046876, p-value = 0.675
## alternative hypothesis: two-sidedshapiro.test(res.est_CCI.25)##
## Shapiro-Wilk normality test
##
## data: res.est_CCI.25
## W = 0.99303, p-value = 0.3316# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_CCI.25, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (CCI.R)^lambda.CCI2_25 ~ N, data = BD2_25)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI2_25 by N
## Bartlett's K-squared = 21.26, df = 3, p-value = 9.296e-05bartlett.test( (CCI.R)^lambda.CCI2_25 ~ P, data = BD2_25)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI2_25 by P
## Bartlett's K-squared = 5.1451, df = 2, p-value = 0.07634# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_CCI.25 ~ BD2_25$N * BD2_25$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 4.0452 2.057e-05 ***
## 225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.25 ~ BD2_25$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 3 7.4792 8.418e-05 ***
## 233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.25 ~ BD2_25$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 3.217 0.04186 *
## 234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.25 ~ BD2_25$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 4.0452 2.057e-05 ***
## 225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test levene - Modelo transformado
leveneTest((BD2_25$CCI.R)^lambda.CCI2_25 ~ BD2_25$N * BD2_25$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 3.1953 0.0004701 ***
## 225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD2_25$CCI.R)^lambda.CCI2_25 ~ BD2_25$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 3 4.7534 0.003087 **
## 233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD2_25$CCI.R)^lambda.CCI2_25 ~ BD2_25$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 2.7746 0.06443 .
## 234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.25 ~ BD2_25$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 3.1953 0.0004701 ***
## 225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD2_25$CCI.R)^lambda.CCI2_25 ~ BD2_25$N)##
## hov: Brown-Forsyth
##
## data: (BD2_25$CCI.R)^lambda.CCI2_25
## F = 4.7534, df:BD2_25$N = 3, df:Residuals = 233, p-value = 0.003087
## alternative hypothesis: variances are not identicalhov((BD2_25$CCI.R)^lambda.CCI2_25 ~ BD2_25$P)##
## hov: Brown-Forsyth
##
## data: (BD2_25$CCI.R)^lambda.CCI2_25
## F = 2.7746, df:BD2_25$P = 2, df:Residuals = 234, p-value = 0.06443
## alternative hypothesis: variances are not identicalhov((BD2_25$CCI.R)^lambda.CCI2_25 ~ BD2_25$TTO)##
## hov: Brown-Forsyth
##
## data: (BD2_25$CCI.R)^lambda.CCI2_25
## F = 3.1953, df:BD2_25$TTO = 11, df:Residuals = 225, p-value = 0.0004701
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 25 DDT)
library(agricolae)
mod2_CCI.25 <- aov( (CCI.R)^lambda.CCI2_25 ~ TTO, data = BD2_25)
par(mfrow=c(1,2))
plot(mod2_CCI.25, which = 1)
plot(mod2_CCI.25, which = 2)
AOV2_CCI.25 <- Anova(mod2_CCI.25, type = "III")
Tuk_CCI.25 <- HSD.test(mod2_CCI.25, "TTO", console=F, alpha = 0.05)
Tukey_BD2.CCI.25 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_CCI.25$groups[2] ), var = "TTO"),
SDT=rep(x=25, times=12))
Tukey_BD2.CCI.25## TTO groups SDT
## 1 1.1N/0.5P a 25
## 2 1.1N/0.0P ab 25
## 3 1.1N/0.25P abc 25
## 4 2.5N/0.25P abcd 25
## 5 2.5N/0.5P bcde 25
## 6 1.8N/0.0P bcde 25
## 7 1.8N/0.5P cde 25
## 8 0.0N/0.0P de 25
## 9 1.8N/0.25P de 25
## 10 2.5N/0.0P de 25
## 11 0.0N/0.5P e 25
## 12 0.0N/0.25P e 25# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD2.CCI.N25 <- kruskal.test(BD2_25$CCI.R, BD2_25$N); KW.BD2.CCI.N25##
## Kruskal-Wallis rank sum test
##
## data: BD2_25$CCI.R and BD2_25$N
## Kruskal-Wallis chi-squared = 77.869, df = 3, p-value < 2.2e-16KW.BD2.CCI.P25 <- kruskal.test(BD2_25$CCI.R, BD2_25$P); KW.BD2.CCI.P25 ##
## Kruskal-Wallis rank sum test
##
## data: BD2_25$CCI.R and BD2_25$P
## Kruskal-Wallis chi-squared = 4.602, df = 2, p-value = 0.1002KW.BD2.CCI.TTO25 <- kruskal.test(BD2_25$CCI.R, BD2_25$TTO); KW.BD2.CCI.TTO25##
## Kruskal-Wallis rank sum test
##
## data: BD2_25$CCI.R and BD2_25$TTO
## Kruskal-Wallis chi-squared = 98.107, df = 11, p-value = 4.23e-16# ver comparaciones puntuales entre tratamientos
KWlt.BD2.CCI.N25 <- kruskal(BD2_25$CCI.R, BD2_25$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.CCI.P25 <- kruskal(BD2_25$CCI.R, BD2_25$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.CCI.TTO25 <- kruskal(BD2_25$CCI.R, BD2_25$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.05)
kruskal_BD2.CCI.25 <- data.frame(
TTO=row.names(KWlt.BD2.CCI.TTO25$groups),
groups_np= c(KWlt.BD2.CCI.TTO25$groups[2])
)
kruskal_BD2.CCI.25## TTO groups
## 1 1.1N/0.5P a
## 2 1.1N/0.0P ab
## 3 2.5N/0.25P abc
## 4 1.1N/0.25P bc
## 5 1.8N/0.0P bc
## 6 2.5N/0.5P bcd
## 7 1.8N/0.5P cd
## 8 0.0N/0.0P cde
## 9 1.8N/0.25P cde
## 10 2.5N/0.0P cde
## 11 0.0N/0.5P de
## 12 0.0N/0.25P e# Grafico CCI 25
Resumen_CCI2_25 <- BD2_25 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(CCI.R), media=mean(CCI.R), sd=sd(CCI.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD2.CCI.25, by="TTO") %>%
dplyr::full_join(., kruskal_BD2.CCI.25, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_CCI2_25, n = 12)## # A tibble: 12 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0.0P 28 16.1 16.1 4.87 0.920 1.89 de 25 cde
## 2 0.0N/0.25P 24 11.8 12.1 4.98 1.02 2.10 e 25 e
## 3 0.0N/0.5P 24 12.8 12.7 4.80 0.980 2.03 e 25 de
## 4 1.1N/0.0P 24 24.7 26.6 7.03 1.44 2.97 ab 25 ab
## 5 1.1N/0.25P 24 23.7 24 11.2 2.29 4.73 abc 25 bc
## 6 1.1N/0.5P 24 27.8 30.5 6.71 1.37 2.83 a 25 a
## 7 1.8N/0.0P 24 22.4 20.0 6.78 1.38 2.86 bcde 25 bc
## 8 1.8N/0.25P 14 16.3 15.6 7.10 1.90 4.10 de 25 cde
## 9 1.8N/0.5P 22 20.5 18.8 9.04 1.93 4.01 cde 25 cd
## 10 2.5N/0.0P 8 10.8 13.0 6.37 2.25 5.33 de 25 cde
## 11 2.5N/0.25P 9 24.8 22.8 6.95 2.32 5.35 abcd 25 abc
## 12 2.5N/0.5P 12 19.2 20.3 7.19 2.08 4.57 bcde 25 bcdFigura_CCI2_25 <- ggplot(Resumen_CCI2_25, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+3.3 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-5.3 ),color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (CCI)")+
ylim(0, 50)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_CCI2_25
### 16 SDT
# Subconjunto base datos
BD2_16 <- as.data.frame(subset(BD2, subset = SDT==16,
select = c(TTO, N, P, CCI.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_CCI2_16 <- data.frame(lam16 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam16 <- seq(from = -10, to = 20, by = 0.5)
for (i in lam16) {
N <- bartlett.test((CCI.R)^i ~ N, data = BD2_16)
P <- bartlett.test((CCI.R)^i ~ P, data = BD2_16)
# Añadir los valores al dataframe
df_CCI2_16 <- rbind(df_CCI2_16, data.frame(lam16 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_CCI2_16$N.p.value~df_CCI2_16$lam16, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_CCI2_16$P.p.value~df_CCI2_16$lam16, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(CCI.R ~ N * P, data = BD2_16, lambda = seq(-10, 20, 1/10))
#lambda.CCI2_16 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.CCI2_16 <- 0.4
# Creación del modelo ANOVA de dos factores tipo III
mod_CCI.16 <- aov( (CCI.R)^lambda.CCI2_16 ~ N * P, data = BD2_16)
AOV_CCI.16 <- Anova(mod_CCI.16, type = "III"); AOV_CCI.16## Anova Table (Type III tests)
##
## Response: (CCI.R)^lambda.CCI2_16
## Sum Sq Df F value Pr(>F)
## (Intercept) 227.423 1 709.0034 < 2.2e-16 ***
## N 3.285 3 3.4142 0.01795 *
## P 2.289 2 3.5681 0.02953 *
## N:P 10.694 6 5.5565 1.877e-05 ***
## Residuals 87.569 273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_CCI.16 <- mod_CCI.16$residuals
res.est_CCI.16 <-res_CCI.16/sd(res_CCI.16)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD2_16$CCI.R), main="Variable original" )
hist( (BD2_16$CCI.R)^lambda.CCI2_16, main="Variable transformada" )
boxplot( (BD2_16$CCI.R), horizontal = T, main="Variable original" )
boxplot( (BD2_16$CCI.R)^lambda.CCI2_16, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_CCI.16, which = 2)
library(nortest)
ad.test(res.est_CCI.16)##
## Anderson-Darling normality test
##
## data: res.est_CCI.16
## A = 1.0454, p-value = 0.009374cvm.test(res.est_CCI.16)##
## Cramer-von Mises normality test
##
## data: res.est_CCI.16
## W = 0.18283, p-value = 0.008826pearson.test(res.est_CCI.16)##
## Pearson chi-square normality test
##
## data: res.est_CCI.16
## P = 42.368, p-value = 0.0005932ks.test(res.est_CCI.16, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_CCI.16
## D = 0.054609, p-value = 0.3632
## alternative hypothesis: two-sidedshapiro.test(res.est_CCI.16)##
## Shapiro-Wilk normality test
##
## data: res.est_CCI.16
## W = 0.98967, p-value = 0.04087# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_CCI.16, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (CCI.R)^lambda.CCI2_16 ~ N, data = BD2_16)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI2_16 by N
## Bartlett's K-squared = 12.146, df = 3, p-value = 0.006898bartlett.test( (CCI.R)^lambda.CCI2_16 ~ P, data = BD2_16)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI2_16 by P
## Bartlett's K-squared = 5.0534, df = 2, p-value = 0.07992# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_CCI.16 ~ BD2_16$N * BD2_16$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 5.7942 1.992e-08 ***
## 273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.16 ~ BD2_16$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 3 3.4399 0.01731 *
## 281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.16 ~ BD2_16$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 6.955 0.001126 **
## 282
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.16 ~ BD2_16$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 5.7942 1.992e-08 ***
## 273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test levene - Modelo transformado
leveneTest((BD2_16$CCI.R)^lambda.CCI2_16 ~ BD2_16$N * BD2_16$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 5.0483 3.551e-07 ***
## 273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD2_16$CCI.R)^lambda.CCI2_16 ~ BD2_16$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 3 2.5526 0.05582 .
## 281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD2_16$CCI.R)^lambda.CCI2_16 ~ BD2_16$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 3.5132 0.03111 *
## 282
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.16 ~ BD2_16$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 5.0483 3.551e-07 ***
## 273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD2_16$CCI.R)^lambda.CCI2_16 ~ BD2_16$N)##
## hov: Brown-Forsyth
##
## data: (BD2_16$CCI.R)^lambda.CCI2_16
## F = 2.5526, df:BD2_16$N = 3, df:Residuals = 281, p-value = 0.05582
## alternative hypothesis: variances are not identicalhov((BD2_16$CCI.R)^lambda.CCI2_16 ~ BD2_16$P)##
## hov: Brown-Forsyth
##
## data: (BD2_16$CCI.R)^lambda.CCI2_16
## F = 3.5132, df:BD2_16$P = 2, df:Residuals = 282, p-value = 0.03111
## alternative hypothesis: variances are not identicalhov((BD2_16$CCI.R)^lambda.CCI2_16 ~ BD2_16$TTO)##
## hov: Brown-Forsyth
##
## data: (BD2_16$CCI.R)^lambda.CCI2_16
## F = 5.0483, df:BD2_16$TTO = 11, df:Residuals = 273, p-value = 3.551e-07
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 16 DDT)
library(agricolae)
mod2_CCI.16 <- aov( (CCI.R)^1 ~ TTO, data = BD2_16)
par(mfrow=c(1,2))
plot(mod2_CCI.16, which = 1)
plot(mod2_CCI.16, which = 2)
AOV2_CCI.16 <- Anova(mod2_CCI.16, type = "III")
Tuk_CCI.16 <- HSD.test(mod2_CCI.16, "TTO", console=F, alpha = 0.05)
Tukey_BD2.CCI.16 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_CCI.16$groups[2] ), var = "TTO"),
SDT=rep(x=16, times=12))
Tukey_BD2.CCI.16## TTO groups SDT
## 1 1.8N/0.5P a 16
## 2 1.1N/0.5P ab 16
## 3 1.1N/0.0P abc 16
## 4 2.5N/0.25P abcd 16
## 5 1.1N/0.25P abcd 16
## 6 1.8N/0.0P abcd 16
## 7 2.5N/0.5P abcde 16
## 8 0.0N/0.25P bcde 16
## 9 2.5N/0.0P cde 16
## 10 0.0N/0.0P cde 16
## 11 1.8N/0.25P de 16
## 12 0.0N/0.5P e 16# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD2.CCI.N16 <- kruskal.test(BD2_16$CCI.R, BD2_16$N); KW.BD2.CCI.N16##
## Kruskal-Wallis rank sum test
##
## data: BD2_16$CCI.R and BD2_16$N
## Kruskal-Wallis chi-squared = 27.242, df = 3, p-value = 5.238e-06KW.BD2.CCI.P16 <- kruskal.test(BD2_16$CCI.R, BD2_16$P); KW.BD2.CCI.P16 ##
## Kruskal-Wallis rank sum test
##
## data: BD2_16$CCI.R and BD2_16$P
## Kruskal-Wallis chi-squared = 2.2782, df = 2, p-value = 0.3201KW.BD2.CCI.TTO16 <- kruskal.test(BD2_16$CCI.R, BD2_16$TTO); KW.BD2.CCI.TTO16##
## Kruskal-Wallis rank sum test
##
## data: BD2_16$CCI.R and BD2_16$TTO
## Kruskal-Wallis chi-squared = 59.934, df = 11, p-value = 9.536e-09# ver comparaciones puntuales entre tratamientos
KWlt.BD2.CCI.N16 <- kruskal(BD2_16$CCI.R, BD2_16$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.CCI.P16 <- kruskal(BD2_16$CCI.R, BD2_16$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.CCI.TTO16 <- kruskal(BD2_16$CCI.R, BD2_16$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.05)
kruskal_BD2.CCI.16 <- data.frame(
TTO=row.names(KWlt.BD2.CCI.TTO16$groups),
groups_np= c(KWlt.BD2.CCI.TTO16$groups[2])
)
kruskal_BD2.CCI.16## TTO groups
## 1 1.1N/0.5P a
## 2 1.8N/0.5P ab
## 3 2.5N/0.25P abc
## 4 1.1N/0.0P abc
## 5 1.8N/0.0P abc
## 6 1.1N/0.25P abc
## 7 2.5N/0.5P abcd
## 8 0.0N/0.25P bcd
## 9 2.5N/0.0P bcd
## 10 0.0N/0.0P cd
## 11 1.8N/0.25P cd
## 12 0.0N/0.5P d# Grafico CCI 16
Resumen_CCI2_16 <- BD2_16 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(CCI.R), media=mean(CCI.R), sd=sd(CCI.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD2.CCI.16, by="TTO") %>%
dplyr::full_join(., kruskal_BD2.CCI.16, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_CCI2_16, n = 12)## # A tibble: 12 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0.0P 24 18.6 17.7 7.85 1.60 3.31 cde 16 cd
## 2 0.0N/0.25P 24 19.7 19.4 3.92 0.800 1.65 bcde 16 bcd
## 3 0.0N/0.5P 24 14.2 13.8 4.13 0.843 1.74 e 16 d
## 4 1.1N/0.0P 24 21.4 24.8 13.4 2.74 5.68 abc 16 abc
## 5 1.1N/0.25P 24 23.1 22.6 8.64 1.76 3.65 abcd 16 abc
## 6 1.1N/0.5P 24 27.7 26.7 3.63 0.741 1.53 ab 16 a
## 7 1.8N/0.0P 24 21.4 22.3 7.42 1.51 3.13 abcd 16 abc
## 8 1.8N/0.25P 24 16.4 16.3 6.88 1.40 2.91 de 16 cd
## 9 1.8N/0.5P 23 24.9 29.0 14.0 2.92 6.05 a 16 ab
## 10 2.5N/0.0P 24 19.6 18.0 8.50 1.73 3.59 cde 16 bcd
## 11 2.5N/0.25P 23 22.4 23.3 9.44 1.97 4.08 abcd 16 abc
## 12 2.5N/0.5P 23 20.2 21.4 9.95 2.07 4.30 abcde 16 abcdFigura_CCI2_16 <- ggplot(Resumen_CCI2_16, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+3.3 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-3.3 ), color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (CCI)")+
ylim(0.0, 60.0)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_CCI2_16
07 SDT
# Subconjunto base datos
BD2_07 <- as.data.frame(subset(BD2, subset = SDT==07,
select = c(TTO, N, P, CCI.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_CCI2_07 <- data.frame(lam07 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam07 <- seq(from = -10, to = 20, by = 0.5)
for (i in lam07) {
N <- bartlett.test((CCI.R)^i ~ N, data = BD2_07)
P <- bartlett.test((CCI.R)^i ~ P, data = BD2_07)
# Añadir los valores al dataframe
df_CCI2_07 <- rbind(df_CCI2_07, data.frame(lam07 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_CCI2_07$N.p.value~df_CCI2_07$lam07, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_CCI2_07$P.p.value~df_CCI2_07$lam07, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(CCI.R ~ N * P, data = BD2_07, lambda = seq(-10, 20, 1/10))
#lambda.CCI2_07 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.CCI2_07 <- 1
# Creación del modelo ANOVA de dos factores tipo III
mod_CCI.07 <- aov( (CCI.R)^lambda.CCI2_07 ~ N * P, data = BD2_07)
AOV_CCI.07 <- Anova(mod_CCI.07, type = "III"); AOV_CCI.07## Anova Table (Type III tests)
##
## Response: (CCI.R)^lambda.CCI2_07
## Sum Sq Df F value Pr(>F)
## (Intercept) 4865.5 1 98.5702 < 2.2e-16 ***
## N 1099.4 3 7.4243 8.512e-05 ***
## P 409.5 2 4.1482 0.01681 *
## N:P 598.8 6 2.0220 0.06298 .
## Residuals 13376.9 271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_CCI.07 <- mod_CCI.07$residuals
res.est_CCI.07 <-res_CCI.07/sd(res_CCI.07)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD2_07$CCI.R), main="Variable original" )
hist( (BD2_07$CCI.R)^lambda.CCI2_07, main="Variable transformada" )
boxplot( (BD2_07$CCI.R), horizontal = T, main="Variable original" )
boxplot( (BD2_07$CCI.R)^lambda.CCI2_07, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_CCI.07, which = 2)
library(nortest)
ad.test(res.est_CCI.07)##
## Anderson-Darling normality test
##
## data: res.est_CCI.07
## A = 1.4314, p-value = 0.001047cvm.test(res.est_CCI.07)##
## Cramer-von Mises normality test
##
## data: res.est_CCI.07
## W = 0.19841, p-value = 0.005584pearson.test(res.est_CCI.07)##
## Pearson chi-square normality test
##
## data: res.est_CCI.07
## P = 25.622, p-value = 0.08163ks.test(res.est_CCI.07, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_CCI.07
## D = 0.056626, p-value = 0.3243
## alternative hypothesis: two-sidedshapiro.test(res.est_CCI.07)##
## Shapiro-Wilk normality test
##
## data: res.est_CCI.07
## W = 0.97328, p-value = 3.908e-05# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_CCI.07, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (CCI.R)^lambda.CCI2_07 ~ N, data = BD2_07)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI2_07 by N
## Bartlett's K-squared = 44.986, df = 3, p-value = 9.315e-10bartlett.test( (CCI.R)^lambda.CCI2_07 ~ P, data = BD2_07)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI2_07 by P
## Bartlett's K-squared = 0.0083638, df = 2, p-value = 0.9958# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_CCI.07 ~ BD2_07$N * BD2_07$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 8.9272 1.479e-13 ***
## 271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.07 ~ BD2_07$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 3 19.124 2.628e-11 ***
## 279
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.07 ~ BD2_07$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 0.4933 0.6111
## 280leveneTest(res.est_CCI.07 ~ BD2_07$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 8.9272 1.479e-13 ***
## 271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test levene - Modelo transformado
leveneTest((BD2_07$CCI.R)^lambda.CCI2_07 ~ BD2_07$N * BD2_07$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 7.4754 3.309e-11 ***
## 271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD2_07$CCI.R)^lambda.CCI2_07 ~ BD2_07$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 3 10.367 1.726e-06 ***
## 279
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD2_07$CCI.R)^lambda.CCI2_07 ~ BD2_07$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 0.0928 0.9114
## 280leveneTest(res.est_CCI.07 ~ BD2_07$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 7.4754 3.309e-11 ***
## 271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD2_07$CCI.R)^lambda.CCI2_07 ~ BD2_07$N)##
## hov: Brown-Forsyth
##
## data: (BD2_07$CCI.R)^lambda.CCI2_07
## F = 10.367, df:BD2_07$N = 3, df:Residuals = 279, p-value = 1.726e-06
## alternative hypothesis: variances are not identicalhov((BD2_07$CCI.R)^lambda.CCI2_07 ~ BD2_07$P)##
## hov: Brown-Forsyth
##
## data: (BD2_07$CCI.R)^lambda.CCI2_07
## F = 0.092807, df:BD2_07$P = 2, df:Residuals = 280, p-value = 0.9114
## alternative hypothesis: variances are not identicalhov((BD2_07$CCI.R)^lambda.CCI2_07 ~ BD2_07$TTO)##
## hov: Brown-Forsyth
##
## data: (BD2_07$CCI.R)^lambda.CCI2_07
## F = 7.4754, df:BD2_07$TTO = 11, df:Residuals = 271, p-value = 3.309e-11
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 07 DDT)
library(agricolae)
mod2_CCI.07 <- aov( (CCI.R)^lambda.CCI2_07 ~ TTO, data = BD2_07)
par(mfrow=c(1,2))
plot(mod2_CCI.07, which = 1)
plot(mod2_CCI.07, which = 2)
AOV2_CCI.07 <- Anova(mod2_CCI.07, type = "III")
Tuk_CCI.07 <- HSD.test(mod2_CCI.07, "TTO", console=F, alpha = 0.05)
Tukey_BD2.CCI.07 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_CCI.07$groups[2] ), var = "TTO"),
SDT=rep(x=07, times=12))
Tukey_BD2.CCI.07## TTO groups SDT
## 1 1.1N/0.0P a 7
## 2 1.1N/0.25P a 7
## 3 2.5N/0.25P a 7
## 4 1.8N/0.0P ab 7
## 5 0.0N/0.25P ab 7
## 6 2.5N/0.5P ab 7
## 7 2.5N/0.0P ab 7
## 8 1.8N/0.5P ab 7
## 9 1.8N/0.25P ab 7
## 10 1.1N/0.5P ab 7
## 11 0.0N/0.0P b 7
## 12 0.0N/0.5P b 7# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD2.CCI.N07 <- kruskal.test(BD2_07$CCI.R, BD2_07$N); KW.BD2.CCI.N07##
## Kruskal-Wallis rank sum test
##
## data: BD2_07$CCI.R and BD2_07$N
## Kruskal-Wallis chi-squared = 13.142, df = 3, p-value = 0.004339KW.BD2.CCI.P07 <- kruskal.test(BD2_07$CCI.R, BD2_07$P); KW.BD2.CCI.P07 ##
## Kruskal-Wallis rank sum test
##
## data: BD2_07$CCI.R and BD2_07$P
## Kruskal-Wallis chi-squared = 9.3126, df = 2, p-value = 0.009502KW.BD2.CCI.TTO07 <- kruskal.test(BD2_07$CCI.R, BD2_07$TTO); KW.BD2.CCI.TTO07##
## Kruskal-Wallis rank sum test
##
## data: BD2_07$CCI.R and BD2_07$TTO
## Kruskal-Wallis chi-squared = 33.234, df = 11, p-value = 0.0004822# ver comparaciones puntuales entre tratamientos
KWlt.BD2.CCI.N07 <- kruskal(BD2_07$CCI.R, BD2_07$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.CCI.P07 <- kruskal(BD2_07$CCI.R, BD2_07$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.CCI.TTO07 <- kruskal(BD2_07$CCI.R, BD2_07$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.05)
kruskal_BD2.CCI.07 <- data.frame(
TTO=row.names(KWlt.BD2.CCI.TTO07$groups),
groups_np= c(KWlt.BD2.CCI.TTO07$groups[2])
)
kruskal_BD2.CCI.07## TTO groups
## 1 1.1N/0.0P a
## 2 2.5N/0.25P a
## 3 1.1N/0.25P ab
## 4 0.0N/0.25P ab
## 5 1.8N/0.0P ab
## 6 2.5N/0.0P ab
## 7 1.8N/0.25P ab
## 8 1.8N/0.5P ab
## 9 1.1N/0.5P ab
## 10 2.5N/0.5P ab
## 11 0.0N/0.0P b
## 12 0.0N/0.5P b# Grafico CCI 07
Resumen_CCI2_07 <- BD2_07 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(CCI.R), media=mean(CCI.R), sd=sd(CCI.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD2.CCI.07, by="TTO") %>%
dplyr::full_join(., kruskal_BD2.CCI.07, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_CCI2_07, n = 12)## # A tibble: 12 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0.0P 28 13.2 13.2 3.60 0.680 1.40 b 7 b
## 2 0.0N/0.25P 24 20.2 18.1 4.08 0.833 1.72 ab 7 ab
## 3 0.0N/0.5P 24 12.9 13.0 2.92 0.596 1.23 b 7 b
## 4 1.1N/0.0P 20 20.2 22.8 11.1 2.47 5.17 a 7 a
## 5 1.1N/0.25P 19 19.9 21.5 9.98 2.29 4.81 a 7 ab
## 6 1.1N/0.5P 24 15.4 16.3 5.42 1.11 2.29 ab 7 ab
## 7 1.8N/0.0P 28 17.6 18.2 4.97 0.939 1.93 ab 7 ab
## 8 1.8N/0.25P 24 17.0 17.3 6.77 1.38 2.86 ab 7 ab
## 9 1.8N/0.5P 24 15.8 17.4 7.67 1.57 3.24 ab 7 ab
## 10 2.5N/0.0P 24 16.8 17.4 6.38 1.30 2.69 ab 7 ab
## 11 2.5N/0.25P 20 20.4 20.3 7.66 1.71 3.58 a 7 a
## 12 2.5N/0.5P 24 16.4 18.1 10.5 2.14 4.44 ab 7 abFigura_CCI2_07 <- ggplot(Resumen_CCI2_07, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+3.3 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-3.3 ), color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (CCI)")+
ylim(0.0, 60.0)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_CCI2_07
par(mfrow=c(1,3))
Figura_CCI2_07
Figura_CCI2_16
Figura_CCI2_25
## Para FC ### Exploratorio FC
ggplot(BD2, aes(x=factor(SDT, levels = c("7", "16", "25")), y=FC.R))+
geom_boxplot(alpha=0.2)+
geom_jitter(aes(color=B), width = 0.1, alpha=0.5)+
facet_grid(~TTO) + ylim(c(0, 1))
### 25 SDT
# Subconjunto base datos
BD2_25 <- as.data.frame(subset(BD2, subset = SDT==25,
select = c(TTO, N, P, FC.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_FC2_25 <- data.frame(lam25 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam25 <- seq(from = -1, to = 20, by = 0.5)
for (i in lam25) {
N <- bartlett.test((FC.R)^i ~ N, data = BD2_25)
P <- bartlett.test((FC.R)^i ~ P, data = BD2_25)
# Añadir los valores al dataframe
df_FC2_25 <- rbind(df_FC2_25, data.frame(lam25 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_FC2_25$N.p.value~df_FC2_25$lam25, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_FC2_25$P.p.value~df_FC2_25$lam25, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(FC.R ~ N * P, data = BD2_25, lambda = seq(-1, 20, 1/10))
#lambda.FC2_25 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.FC2_25 <- 8
# Creación del modelo ANOVA de dos factores tipo III
mod_FC.25 <- aov( (FC.R)^lambda.FC2_25 ~ N * P, data = BD2_25)
AOV_FC.25 <- Anova(mod_FC.25, type = "III"); AOV_FC.25## Anova Table (Type III tests)
##
## Response: (FC.R)^lambda.FC2_25
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.39258 1 211.1659 < 2.2e-16 ***
## N 0.06156 3 11.0383 8.662e-07 ***
## P 0.00201 2 0.5415 0.58262
## N:P 0.02333 6 2.0917 0.05519 .
## Residuals 0.41829 225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_FC.25 <- mod_FC.25$residuals
res.est_FC.25 <-res_FC.25/sd(res_FC.25)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD2_25$FC.R), main="Variable original" )
hist( (BD2_25$FC.R)^lambda.FC2_25, main="Variable transformada" )
boxplot( (BD2_25$FC.R), horizontal = T, main="Variable original" )
boxplot( (BD2_25$FC.R)^lambda.FC2_25, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_FC.25, which = 2)
library(nortest)
ad.test(res.est_FC.25)##
## Anderson-Darling normality test
##
## data: res.est_FC.25
## A = 2.6541, p-value = 1.041e-06cvm.test(res.est_FC.25)##
## Cramer-von Mises normality test
##
## data: res.est_FC.25
## W = 0.41947, p-value = 1.645e-05pearson.test(res.est_FC.25)##
## Pearson chi-square normality test
##
## data: res.est_FC.25
## P = 36.797, p-value = 0.001355ks.test(res.est_FC.25, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_FC.25
## D = 0.083481, p-value = 0.07352
## alternative hypothesis: two-sidedshapiro.test(res.est_FC.25)##
## Shapiro-Wilk normality test
##
## data: res.est_FC.25
## W = 0.96171, p-value = 5.715e-06# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_FC.25, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (FC.R)^lambda.FC2_25 ~ N, data = BD2_25)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC2_25 by N
## Bartlett's K-squared = 1.8854, df = 3, p-value = 0.5965bartlett.test( (FC.R)^lambda.FC2_25 ~ P, data = BD2_25)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC2_25 by P
## Bartlett's K-squared = 6.2515, df = 2, p-value = 0.0439# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_FC.25 ~ BD2_25$N * BD2_25$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 3.8998 3.53e-05 ***
## 225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.25 ~ BD2_25$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 3 1.5295 0.2076
## 233leveneTest(res.est_FC.25 ~ BD2_25$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 12.659 6.023e-06 ***
## 234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.25 ~ BD2_25$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 3.8998 3.53e-05 ***
## 225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test levene - Modelo transformado
leveneTest((BD2_25$FC.R)^lambda.FC2_25 ~ BD2_25$N * BD2_25$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 2.1827 0.01622 *
## 225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD2_25$FC.R)^lambda.FC2_25 ~ BD2_25$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 3 0.6198 0.6028
## 233leveneTest((BD2_25$FC.R)^lambda.FC2_25 ~ BD2_25$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 3.0543 0.04904 *
## 234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.25 ~ BD2_25$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 2.1827 0.01622 *
## 225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD2_25$FC.R)^lambda.FC2_25 ~ BD2_25$N)##
## hov: Brown-Forsyth
##
## data: (BD2_25$FC.R)^lambda.FC2_25
## F = 0.61985, df:BD2_25$N = 3, df:Residuals = 233, p-value = 0.6028
## alternative hypothesis: variances are not identicalhov((BD2_25$FC.R)^lambda.FC2_25 ~ BD2_25$P)##
## hov: Brown-Forsyth
##
## data: (BD2_25$FC.R)^lambda.FC2_25
## F = 3.0543, df:BD2_25$P = 2, df:Residuals = 234, p-value = 0.04904
## alternative hypothesis: variances are not identicalhov((BD2_25$FC.R)^lambda.FC2_25 ~ BD2_25$TTO)##
## hov: Brown-Forsyth
##
## data: (BD2_25$FC.R)^lambda.FC2_25
## F = 2.1827, df:BD2_25$TTO = 11, df:Residuals = 225, p-value = 0.01622
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 25 DDT)
library(agricolae)
mod2_FC.25 <- aov( (FC.R)^lambda.FC2_25 ~ TTO, data = BD2_25)
par(mfrow=c(1,2))
plot(mod2_FC.25, which = 1)
plot(mod2_FC.25, which = 2)
AOV2_FC.25 <- Anova(mod2_FC.25, type = "III")
Tuk_FC.25 <- HSD.test(mod2_FC.25, "TTO", console=F, alpha = 0.05)
Tukey_BD2.FC.25 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_FC.25$groups[2] ), var = "TTO"),
SDT=rep(x=25, times=12))
Tukey_BD2.FC.25## TTO groups SDT
## 1 0.0N/0.5P a 25
## 2 0.0N/0.25P ab 25
## 3 0.0N/0.0P ab 25
## 4 2.5N/0.5P abc 25
## 5 1.1N/0.0P bcd 25
## 6 1.1N/0.5P bcd 25
## 7 1.1N/0.25P cd 25
## 8 1.8N/0.0P cd 25
## 9 2.5N/0.25P cd 25
## 10 1.8N/0.5P d 25
## 11 1.8N/0.25P d 25
## 12 2.5N/0.0P d 25# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD2.FC.N25 <- kruskal.test(BD2_25$FC.R, BD2_25$N); KW.BD2.FC.N25##
## Kruskal-Wallis rank sum test
##
## data: BD2_25$FC.R and BD2_25$N
## Kruskal-Wallis chi-squared = 68.341, df = 3, p-value = 9.671e-15KW.BD2.FC.P25 <- kruskal.test(BD2_25$FC.R, BD2_25$P); KW.BD2.FC.P25 ##
## Kruskal-Wallis rank sum test
##
## data: BD2_25$FC.R and BD2_25$P
## Kruskal-Wallis chi-squared = 1.2902, df = 2, p-value = 0.5246KW.BD2.FC.TTO25 <- kruskal.test(BD2_25$FC.R, BD2_25$TTO); KW.BD2.FC.TTO25##
## Kruskal-Wallis rank sum test
##
## data: BD2_25$FC.R and BD2_25$TTO
## Kruskal-Wallis chi-squared = 80.791, df = 11, p-value = 1.038e-12# ver comparaciones puntuales entre tratamientos
KWlt.BD2.FC.N25 <- kruskal(BD2_25$FC.R, BD2_25$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.FC.P25 <- kruskal(BD2_25$FC.R, BD2_25$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.FC.TTO25 <- kruskal(BD2_25$FC.R, BD2_25$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.05)
kruskal_BD2.FC.25 <- data.frame(
TTO=row.names(KWlt.BD2.FC.TTO25$groups),
groups_np= c(KWlt.BD2.FC.TTO25$groups[2])
)
kruskal_BD2.FC.25## TTO groups
## 1 0.0N/0.5P a
## 2 0.0N/0.25P a
## 3 0.0N/0.0P ab
## 4 2.5N/0.5P abc
## 5 1.1N/0.0P abcd
## 6 1.1N/0.5P bcd
## 7 1.1N/0.25P cd
## 8 1.8N/0.0P cd
## 9 2.5N/0.25P cd
## 10 1.8N/0.5P d
## 11 1.8N/0.25P d
## 12 2.5N/0.0P d# Grafico FC 25
Resumen_FC2_25 <- BD2_25 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(FC.R), media=mean(FC.R), sd=sd(FC.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD2.FC.25, by="TTO") %>%
dplyr::full_join(., kruskal_BD2.FC.25, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_FC2_25, n = 12)## # A tibble: 12 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0… 28 0.760 0.759 0.0401 0.00757 0.0155 ab 25 ab
## 2 0.0N/0… 24 0.763 0.766 0.0212 0.00433 0.00897 ab 25 a
## 3 0.0N/0… 24 0.774 0.771 0.0298 0.00607 0.0126 a 25 a
## 4 1.1N/0… 24 0.748 0.725 0.0555 0.0113 0.0234 bcd 25 abcd
## 5 1.1N/0… 24 0.691 0.699 0.0538 0.0110 0.0227 cd 25 cd
## 6 1.1N/0… 24 0.722 0.712 0.0663 0.0135 0.0280 bcd 25 bcd
## 7 1.8N/0… 24 0.664 0.663 0.0970 0.0198 0.0410 cd 25 cd
## 8 1.8N/0… 14 0.678 0.675 0.0315 0.00843 0.0182 d 25 d
## 9 1.8N/0… 22 0.626 0.644 0.102 0.0218 0.0454 d 25 d
## 10 2.5N/0… 8 0.624 0.632 0.0919 0.0325 0.0768 d 25 d
## 11 2.5N/0… 9 0.714 0.688 0.0593 0.0198 0.0456 cd 25 cd
## 12 2.5N/0… 12 0.749 0.747 0.0401 0.0116 0.0255 abc 25 abcFigura_FC2_25 <- ggplot(Resumen_FC2_25, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+0.06 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-0.06 ), color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (FC)")+
ylim(0.2, 1.0)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_FC2_25
### 16 SDT
# Subconjunto base datos
BD2_16 <- as.data.frame(subset(BD2, subset = SDT==16,
select = c(TTO, N, P, FC.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_FC2_16 <- data.frame(lam16 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam16 <- seq(from = -10, to = 20, by = 0.5)
for (i in lam16) {
N <- bartlett.test((FC.R)^i ~ N, data = BD2_16)
P <- bartlett.test((FC.R)^i ~ P, data = BD2_16)
# Añadir los valores al dataframe
df_FC2_16 <- rbind(df_FC2_16, data.frame(lam16 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_FC2_16$N.p.value~df_FC2_16$lam16, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_FC2_16$P.p.value~df_FC2_16$lam16, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(FC.R ~ N * P, data = BD2_16, lambda = seq(-10, 20, 1/10))
#lambda.FC2_16 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.FC2_16 <- 7
# Creación del modelo ANOVA de dos factores tipo III
mod_FC.16 <- aov( (FC.R)^lambda.FC2_16 ~ N * P, data = BD2_16)
AOV_FC.16 <- Anova(mod_FC.16, type = "III"); AOV_FC.16## Anova Table (Type III tests)
##
## Response: (FC.R)^lambda.FC2_16
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.8136 1 57.7833 4.717e-13 ***
## N 0.0177 3 0.4193 0.7393
## P 0.0172 2 0.6108 0.5436
## N:P 0.1104 6 1.3065 0.2543
## Residuals 3.8439 273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_FC.16 <- mod_FC.16$residuals
res.est_FC.16 <-res_FC.16/sd(res_FC.16)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD2_16$FC.R), main="Variable original" )
hist( (BD2_16$FC.R)^lambda.FC2_16, main="Variable transformada" )
boxplot( (BD2_16$FC.R), horizontal = T, main="Variable original" )
boxplot( (BD2_16$FC.R)^lambda.FC2_16, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_FC.16, which = 2)
library(nortest)
ad.test(res.est_FC.16)##
## Anderson-Darling normality test
##
## data: res.est_FC.16
## A = 7.0558, p-value < 2.2e-16cvm.test(res.est_FC.16)##
## Cramer-von Mises normality test
##
## data: res.est_FC.16
## W = 1.118, p-value = 7.37e-10pearson.test(res.est_FC.16)##
## Pearson chi-square normality test
##
## data: res.est_FC.16
## P = 87.14, p-value = 2.014e-11ks.test(res.est_FC.16, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_FC.16
## D = 0.12445, p-value = 0.0002933
## alternative hypothesis: two-sidedshapiro.test(res.est_FC.16)##
## Shapiro-Wilk normality test
##
## data: res.est_FC.16
## W = 0.92185, p-value = 4.619e-11# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_FC.16, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (FC.R)^lambda.FC2_16 ~ N, data = BD2_16)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC2_16 by N
## Bartlett's K-squared = 12.016, df = 3, p-value = 0.00733bartlett.test( (FC.R)^lambda.FC2_16 ~ P, data = BD2_16)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC2_16 by P
## Bartlett's K-squared = 0.0076723, df = 2, p-value = 0.9962# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_FC.16 ~ BD2_16$N * BD2_16$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 4.615 1.893e-06 ***
## 273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.16 ~ BD2_16$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 3 6.4892 0.0002928 ***
## 281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.16 ~ BD2_16$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 0.1072 0.8983
## 282leveneTest(res.est_FC.16 ~ BD2_16$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 4.615 1.893e-06 ***
## 273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test levene - Modelo transformado
leveneTest((BD2_16$FC.R)^lambda.FC2_16 ~ BD2_16$N * BD2_16$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 3.108 0.0005808 ***
## 273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD2_16$FC.R)^lambda.FC2_16 ~ BD2_16$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 3 4.5062 0.004171 **
## 281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD2_16$FC.R)^lambda.FC2_16 ~ BD2_16$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 0.0365 0.9642
## 282leveneTest(res.est_FC.16 ~ BD2_16$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 3.108 0.0005808 ***
## 273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD2_16$FC.R)^lambda.FC2_16 ~ BD2_16$N)##
## hov: Brown-Forsyth
##
## data: (BD2_16$FC.R)^lambda.FC2_16
## F = 4.5062, df:BD2_16$N = 3, df:Residuals = 281, p-value = 0.004171
## alternative hypothesis: variances are not identicalhov((BD2_16$FC.R)^lambda.FC2_16 ~ BD2_16$P)##
## hov: Brown-Forsyth
##
## data: (BD2_16$FC.R)^lambda.FC2_16
## F = 0.036476, df:BD2_16$P = 2, df:Residuals = 282, p-value = 0.9642
## alternative hypothesis: variances are not identicalhov((BD2_16$FC.R)^lambda.FC2_16 ~ BD2_16$TTO)##
## hov: Brown-Forsyth
##
## data: (BD2_16$FC.R)^lambda.FC2_16
## F = 3.108, df:BD2_16$TTO = 11, df:Residuals = 273, p-value = 0.0005808
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 16 DDT)
library(agricolae)
mod2_FC.16 <- aov( (FC.R)^lambda.FC2_16 ~ TTO, data = BD2_16)
par(mfrow=c(1,2))
plot(mod2_FC.16, which = 1)
plot(mod2_FC.16, which = 2)
AOV2_FC.16 <- Anova(mod2_FC.16, type = "III")
Tuk_FC.16 <- HSD.test(mod2_FC.16, "TTO", console=F, alpha = 0.05)
Tukey_BD2.FC.16 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_FC.16$groups[2] ), var = "TTO"),
SDT=rep(x=16, times=12))
Tukey_BD2.FC.16## TTO groups SDT
## 1 1.1N/0.25P a 16
## 2 0.0N/0.5P a 16
## 3 2.5N/0.5P a 16
## 4 0.0N/0.25P a 16
## 5 0.0N/0.0P a 16
## 6 1.1N/0.5P a 16
## 7 1.8N/0.5P a 16
## 8 1.1N/0.0P a 16
## 9 1.8N/0.0P a 16
## 10 2.5N/0.0P a 16
## 11 2.5N/0.25P a 16
## 12 1.8N/0.25P a 16# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD2.FC.N16 <- kruskal.test(BD2_16$FC.R, BD2_16$N); KW.BD2.FC.N16##
## Kruskal-Wallis rank sum test
##
## data: BD2_16$FC.R and BD2_16$N
## Kruskal-Wallis chi-squared = 14.177, df = 3, p-value = 0.002674KW.BD2.FC.P16 <- kruskal.test(BD2_16$FC.R, BD2_16$P); KW.BD2.FC.P16 ##
## Kruskal-Wallis rank sum test
##
## data: BD2_16$FC.R and BD2_16$P
## Kruskal-Wallis chi-squared = 5.1452, df = 2, p-value = 0.07634KW.BD2.FC.TTO16 <- kruskal.test(BD2_16$FC.R, BD2_16$TTO); KW.BD2.FC.TTO16##
## Kruskal-Wallis rank sum test
##
## data: BD2_16$FC.R and BD2_16$TTO
## Kruskal-Wallis chi-squared = 24.142, df = 11, p-value = 0.01215# ver comparaciones puntuales entre tratamientos
KWlt.BD2.FC.N16 <- kruskal(BD2_16$FC.R, BD2_16$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.FC.P16 <- kruskal(BD2_16$FC.R, BD2_16$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.FC.TTO16 <- kruskal(BD2_16$FC.R, BD2_16$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.10)
kruskal_BD2.FC.16 <- data.frame(
TTO=row.names(KWlt.BD2.FC.TTO16$groups),
groups_np= c(KWlt.BD2.FC.TTO16$groups[2])
)
kruskal_BD2.FC.16## TTO groups
## 1 0.0N/0.5P a
## 2 1.1N/0.25P ab
## 3 2.5N/0.5P ab
## 4 1.1N/0.5P ab
## 5 0.0N/0.25P ab
## 6 1.1N/0.0P ab
## 7 0.0N/0.0P ab
## 8 1.8N/0.5P ab
## 9 1.8N/0.0P ab
## 10 2.5N/0.0P ab
## 11 2.5N/0.25P ab
## 12 1.8N/0.25P b# Grafico FC 16
Resumen_FC2_16 <- BD2_16 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(FC.R), media=mean(FC.R), sd=sd(FC.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD2.FC.16, by="TTO") %>%
dplyr::full_join(., kruskal_BD2.FC.16, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_FC2_16, n = 12)## # A tibble: 12 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0.… 24 0.762 0.766 0.0711 0.0145 0.0300 a 16 ab
## 2 0.0N/0.… 24 0.768 0.780 0.0480 0.00980 0.0203 a 16 ab
## 3 0.0N/0.… 24 0.778 0.797 0.0466 0.00952 0.0197 a 16 a
## 4 1.1N/0.… 24 0.776 0.772 0.0354 0.00724 0.0150 a 16 ab
## 5 1.1N/0.… 24 0.788 0.784 0.0838 0.0171 0.0354 a 16 ab
## 6 1.1N/0.… 24 0.775 0.778 0.0378 0.00772 0.0160 a 16 ab
## 7 1.8N/0.… 24 0.740 0.722 0.126 0.0257 0.0531 a 16 ab
## 8 1.8N/0.… 24 0.714 0.698 0.123 0.0250 0.0518 a 16 b
## 9 1.8N/0.… 23 0.751 0.728 0.133 0.0277 0.0574 a 16 ab
## 10 2.5N/0.… 24 0.738 0.703 0.133 0.0272 0.0564 a 16 ab
## 11 2.5N/0.… 23 0.747 0.721 0.0917 0.0191 0.0397 a 16 ab
## 12 2.5N/0.… 23 0.775 0.780 0.0805 0.0168 0.0348 a 16 abFigura_FC2_16 <- ggplot(Resumen_FC2_16, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+0.06 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-0.06 ), color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (FC)")+
ylim(0.2, 1.0)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_FC2_16
07 SDT
# Subconjunto base datos
BD2_07 <- as.data.frame(subset(BD2, subset = SDT==07,
select = c(TTO, N, P, FC.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_FC2_07 <- data.frame(lam07 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam07 <- seq(from = -10, to = 20, by = 0.5)
for (i in lam07) {
N <- bartlett.test((FC.R)^i ~ N, data = BD2_07)
P <- bartlett.test((FC.R)^i ~ P, data = BD2_07)
# Añadir los valores al dataframe
df_FC2_07 <- rbind(df_FC2_07, data.frame(lam07 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_FC2_07$N.p.value~df_FC2_07$lam07, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_FC2_07$P.p.value~df_FC2_07$lam07, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(FC.R ~ N * P, data = BD2_07, lambda = seq(-10, 20, 1/10))
#lambda.FC2_07 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.FC2_07 <- 2.8
# Creación del modelo ANOVA de dos factores tipo III
mod_FC.07 <- aov( (FC.R)^lambda.FC2_07 ~ N * P, data = BD2_07)
AOV_FC.07 <- Anova(mod_FC.07, type = "III"); AOV_FC.07## Anova Table (Type III tests)
##
## Response: (FC.R)^lambda.FC2_07
## Sum Sq Df F value Pr(>F)
## (Intercept) 8.7064 1 958.1795 < 2.2e-16 ***
## N 0.2661 3 9.7602 3.896e-06 ***
## P 0.1408 2 7.7480 0.0005343 ***
## N:P 0.3904 6 7.1601 4.223e-07 ***
## Residuals 2.4624 271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_FC.07 <- mod_FC.07$residuals
res.est_FC.07 <-res_FC.07/sd(res_FC.07)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD2_07$FC.R), main="Variable original" )
hist( (BD2_07$FC.R)^lambda.FC2_07, main="Variable transformada" )
boxplot( (BD2_07$FC.R), horizontal = T, main="Variable original" )
boxplot( (BD2_07$FC.R)^lambda.FC2_07, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_FC.07, which = 2)
library(nortest)
ad.test(res.est_FC.07)##
## Anderson-Darling normality test
##
## data: res.est_FC.07
## A = 3.9856, p-value = 6.061e-10cvm.test(res.est_FC.07)##
## Cramer-von Mises normality test
##
## data: res.est_FC.07
## W = 0.78349, p-value = 1.752e-08pearson.test(res.est_FC.07)##
## Pearson chi-square normality test
##
## data: res.est_FC.07
## P = 84.703, p-value = 5.539e-11ks.test(res.est_FC.07, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_FC.07
## D = 0.11502, p-value = 0.001119
## alternative hypothesis: two-sidedshapiro.test(res.est_FC.07)##
## Shapiro-Wilk normality test
##
## data: res.est_FC.07
## W = 0.96731, p-value = 4.891e-06# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_FC.07, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (FC.R)^lambda.FC2_07 ~ N, data = BD2_07)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC2_07 by N
## Bartlett's K-squared = 11.494, df = 3, p-value = 0.009334bartlett.test( (FC.R)^lambda.FC2_07 ~ P, data = BD2_07)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC2_07 by P
## Bartlett's K-squared = 2.2601, df = 2, p-value = 0.323# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_FC.07 ~ BD2_07$N * BD2_07$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 5.2029 1.978e-07 ***
## 271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.07 ~ BD2_07$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 3 5.3446 0.001358 **
## 279
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.07 ~ BD2_07$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 3.3328 0.03711 *
## 280
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.07 ~ BD2_07$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 11 5.2029 1.978e-07 ***
## 271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test levene - Modelo transformado
leveneTest((BD2_07$FC.R)^lambda.FC2_07 ~ BD2_07$N * BD2_07$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 4.598 2.039e-06 ***
## 271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD2_07$FC.R)^lambda.FC2_07 ~ BD2_07$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 3 4.4164 0.004708 **
## 279
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD2_07$FC.R)^lambda.FC2_07 ~ BD2_07$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 0.2535 0.7763
## 280leveneTest(res.est_FC.07 ~ BD2_07$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 11 4.598 2.039e-06 ***
## 271
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD2_07$FC.R)^lambda.FC2_07 ~ BD2_07$N)##
## hov: Brown-Forsyth
##
## data: (BD2_07$FC.R)^lambda.FC2_07
## F = 4.4164, df:BD2_07$N = 3, df:Residuals = 279, p-value = 0.004708
## alternative hypothesis: variances are not identicalhov((BD2_07$FC.R)^lambda.FC2_07 ~ BD2_07$P)##
## hov: Brown-Forsyth
##
## data: (BD2_07$FC.R)^lambda.FC2_07
## F = 0.25348, df:BD2_07$P = 2, df:Residuals = 280, p-value = 0.7763
## alternative hypothesis: variances are not identicalhov((BD2_07$FC.R)^lambda.FC2_07 ~ BD2_07$TTO)##
## hov: Brown-Forsyth
##
## data: (BD2_07$FC.R)^lambda.FC2_07
## F = 4.598, df:BD2_07$TTO = 11, df:Residuals = 271, p-value = 2.039e-06
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 07 DDT)
library(agricolae)
mod2_FC.07 <- aov( (FC.R)^lambda.FC2_07 ~ TTO, data = BD2_07)
par(mfrow=c(1,2))
plot(mod2_FC.07, which = 1)
plot(mod2_FC.07, which = 2)
AOV2_FC.07 <- Anova(mod2_FC.07, type = "III")
Tuk_FC.07 <- HSD.test(mod2_FC.07, "TTO", console=F, alpha = 0.05)
Tukey_BD2.FC.07 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_FC.07$groups[2] ), var = "TTO"),
SDT=rep(x=07, times=12))
Tukey_BD2.FC.07## TTO groups SDT
## 1 0.0N/0.0P a 7
## 2 1.1N/0.25P ab 7
## 3 1.8N/0.5P ab 7
## 4 1.8N/0.0P ab 7
## 5 0.0N/0.5P ab 7
## 6 2.5N/0.25P abc 7
## 7 1.8N/0.25P abc 7
## 8 2.5N/0.0P abc 7
## 9 1.1N/0.5P bc 7
## 10 0.0N/0.25P bc 7
## 11 2.5N/0.5P c 7
## 12 1.1N/0.0P c 7# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD2.FC.N07 <- kruskal.test(BD2_07$FC.R, BD2_07$N); KW.BD2.FC.N07##
## Kruskal-Wallis rank sum test
##
## data: BD2_07$FC.R and BD2_07$N
## Kruskal-Wallis chi-squared = 20.533, df = 3, p-value = 0.0001316KW.BD2.FC.P07 <- kruskal.test(BD2_07$FC.R, BD2_07$P); KW.BD2.FC.P07 ##
## Kruskal-Wallis rank sum test
##
## data: BD2_07$FC.R and BD2_07$P
## Kruskal-Wallis chi-squared = 0.63894, df = 2, p-value = 0.7265KW.BD2.FC.TTO07 <- kruskal.test(BD2_07$FC.R, BD2_07$TTO); KW.BD2.FC.TTO07##
## Kruskal-Wallis rank sum test
##
## data: BD2_07$FC.R and BD2_07$TTO
## Kruskal-Wallis chi-squared = 64.668, df = 11, p-value = 1.244e-09# ver comparaciones puntuales entre tratamientos
KWlt.BD2.FC.N07 <- kruskal(BD2_07$FC.R, BD2_07$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.FC.P07 <- kruskal(BD2_07$FC.R, BD2_07$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD2.FC.TTO07 <- kruskal(BD2_07$FC.R, BD2_07$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.05)
kruskal_BD2.FC.07 <- data.frame(
TTO=row.names(KWlt.BD2.FC.TTO07$groups),
groups_np= c(KWlt.BD2.FC.TTO07$groups[2])
)
kruskal_BD2.FC.07## TTO groups
## 1 0.0N/0.0P a
## 2 0.0N/0.5P ab
## 3 1.8N/0.5P ab
## 4 1.1N/0.25P ab
## 5 1.8N/0.0P ab
## 6 2.5N/0.25P abc
## 7 1.8N/0.25P abc
## 8 2.5N/0.0P abc
## 9 1.1N/0.5P bc
## 10 0.0N/0.25P bc
## 11 2.5N/0.5P c
## 12 1.1N/0.0P c# Grafico FC 07
Resumen_FC2_07 <- BD2_07 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(FC.R), media=mean(FC.R), sd=sd(FC.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD2.FC.07, by="TTO") %>%
dplyr::full_join(., kruskal_BD2.FC.07, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_FC2_07, n = 12)## # A tibble: 12 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0… 28 0.796 0.809 0.0513 0.00970 0.0199 a 7 a
## 2 0.0N/0… 24 0.759 0.754 0.0332 0.00677 0.0140 bc 7 bc
## 3 0.0N/0… 24 0.797 0.794 0.0394 0.00804 0.0166 ab 7 ab
## 4 1.1N/0… 20 0.72 0.732 0.0299 0.00669 0.0140 c 7 c
## 5 1.1N/0… 19 0.789 0.798 0.0655 0.0150 0.0316 ab 7 ab
## 6 1.1N/0… 24 0.756 0.757 0.0700 0.0143 0.0296 bc 7 bc
## 7 1.8N/0… 28 0.793 0.795 0.0605 0.0114 0.0235 ab 7 ab
## 8 1.8N/0… 24 0.776 0.764 0.0746 0.0152 0.0315 abc 7 abc
## 9 1.8N/0… 24 0.801 0.797 0.0575 0.0117 0.0243 ab 7 ab
## 10 2.5N/0… 24 0.763 0.766 0.0228 0.00466 0.00965 abc 7 abc
## 11 2.5N/0… 20 0.767 0.777 0.0703 0.0157 0.0329 abc 7 abc
## 12 2.5N/0… 24 0.739 0.737 0.0254 0.00519 0.0107 c 7 cFigura_FC2_07 <- ggplot(Resumen_FC2_07, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+0.06 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-0.06 ), color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (FC)")+
ylim(0.2, 1.0)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_FC2_07
par(mfrow=c(1,3))
Figura_FC2_07
Figura_FC2_16
Figura_FC2_25
## M D A 2 ### Exploratorio MDA
MDA_Agraz2 <- read_excel("MDA Agraz.xlsx", sheet = "MDA ENSAYO 1")
head(MDA_Agraz2)## # A tibble: 6 × 21
## SDT TTO N P PL M V `440-` `532-` `600-` `440+` `532+`
## <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 25 1 0.00 N 0.00… 101 0.312 5 0.521 0.507 0.275 1.00 0.758
## 2 25 1 0.00 N 0.00… 104 0.315 5 0.718 0.58 0.234 1.65 1.36
## 3 25 1 0.00 N 0.00… 207 0.309 5 0.438 0.363 0.11 1.10 0.934
## 4 25 2 1.10 N 0.00… 9 0.305 5 0.168 0.135 0.124 0.666 0.367
## 5 25 2 1.10 N 0.00… NA NA NA NA NA NA NA NA
## 6 25 2 1.10 N 0.00… 13 0.318 5 0.246 0.261 0.133 0.616 0.403
## # ℹ 9 more variables: `600+` <dbl>, A <dbl>, B <dbl>, C <dbl>, M1 <dbl>,
## # M1.1 <dbl>, M2 <dbl>, M2.1 <dbl>, M2.2 <dbl>MDA_Agraz2 <- MDA_Agraz2 %>%
unite(., col="TTO", c("N","P"), sep = "/", remove = FALSE)
Resumen2 <- MDA_Agraz2 %>%
dplyr::group_by(., N) %>%
dplyr::summarise(.,
mediaM2=mean(M2.2), medianaM2=median(M2.2), stdM2=sd(M2.2), n2= n(),
mediaM1=mean(M1.1), medianaM1=median(M1.1), stdM1=sd(M1.1), n1= n() )
Resumen2## # A tibble: 4 × 9
## N mediaM2 medianaM2 stdM2 n2 mediaM1 medianaM1 stdM1 n1
## <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <int>
## 1 0.00 N 67.8 69.9 24.0 9 30.0 27.5 15.2 9
## 2 1.10 N 15.4 13 8.86 9 6.92 4.10 8.35 9
## 3 1.80 N 13.8 13 4.74 9 0.756 0.5 4.83 9
## 4 2.50 N 13.7 14 3.99 9 6.24 6 4.21 9ggplot(MDA_Agraz2, aes(x=factor(TTO), y=M2.2) )+
geom_point()+
geom_col(alpha=0.2)+
coord_flip()
ggplot(MDA_Agraz2, aes(x=factor(N), y=M2.2) )+
geom_point()+
geom_col(alpha=0.2)+
coord_flip()
ggplot(MDA_Agraz2, aes(x=factor(P), y=M2.2) )+
geom_point()+
geom_col(alpha=0.2)+
coord_flip()
# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_MDA2 <- data.frame(lam06 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam2 <- seq(from = -10, to = 20, by = 0.5)
for (i in lam2) {
N <- bartlett.test((M2.2)^i ~ N, data = MDA_Agraz2)
P <- bartlett.test((M2.2)^i ~ P, data = MDA_Agraz2)
# Añadir los valores al dataframe
df_MDA2 <- rbind(df_MDA2, data.frame(lam2 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_MDA2$N.p.value~df_MDA2$lam2, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_MDA2$P.p.value~df_MDA2$lam2, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
# b <- boxcox(M2.2 ~ N * P, data = MDA_Agraz2, lambda = seq(-10, 20, 1/10))
# lambda.MDA2 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.MDA2 <- 0.1
### Anova MDA
# Creación del modelo ANOVA de dos factores tipo III
mod_MDA2 <- aov( (M2.2)^lambda.MDA2 ~ N * P, data = MDA_Agraz2)
AOV_MDA2 <- Anova(mod_MDA2, type = "III"); AOV_MDA2## Anova Table (Type III tests)
##
## Response: (M2.2)^lambda.MDA2
## Sum Sq Df F value Pr(>F)
## (Intercept) 6.7628 1 2406.8333 < 2.2e-16 ***
## N 0.1036 3 12.2870 4.543e-05 ***
## P 0.0021 2 0.3730 0.6926
## N:P 0.0082 6 0.4887 0.8102
## Residuals 0.0674 24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_MDA2 <- mod_MDA2$residuals
res.est_MDA2 <-res_MDA2/sd(res_MDA2)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (MDA_Agraz2$M2.2), main="Variable original" )
hist( (MDA_Agraz2$M2.2)^lambda.MDA2, main="Variable transformada" )
boxplot( (MDA_Agraz2$M2.2), horizontal = T, main="Variable original" )
boxplot( (MDA_Agraz2$M2.2)^lambda.MDA2, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_MDA2, which = 2)
library(nortest)
ad.test(res.est_MDA2)##
## Anderson-Darling normality test
##
## data: res.est_MDA2
## A = 0.40309, p-value = 0.3394cvm.test(res.est_MDA2)##
## Cramer-von Mises normality test
##
## data: res.est_MDA2
## W = 0.059669, p-value = 0.3727pearson.test(res.est_MDA2)##
## Pearson chi-square normality test
##
## data: res.est_MDA2
## P = 13, p-value = 0.04304ks.test(res.est_MDA2, "pnorm")##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: res.est_MDA2
## D = 0.10972, p-value = 0.7381
## alternative hypothesis: two-sidedshapiro.test(res.est_MDA2)##
## Shapiro-Wilk normality test
##
## data: res.est_MDA2
## W = 0.97236, p-value = 0.4935# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_MDA2, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (M2.2)^lambda.MDA2 ~ N, data = MDA_Agraz2)##
## Bartlett test of homogeneity of variances
##
## data: (M2.2)^lambda.MDA2 by N
## Bartlett's K-squared = 2.5185, df = 3, p-value = 0.472bartlett.test( (M2.2)^lambda.MDA2 ~ P, data = MDA_Agraz2)##
## Bartlett test of homogeneity of variances
##
## data: (M2.2)^lambda.MDA2 by P
## Bartlett's K-squared = 0.43369, df = 2, p-value = 0.8051# Analisis parametrico (Tukey 06 DDT)
library(agricolae)
Tuk_MD2 <- HSD.test(mod_MDA2, "N", console=F, alpha = 0.05)
Tukey_MDA2 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_MD2$groups[2] ), var = "N")
)
Tukey_MDA2$N <- factor(Tukey_MDA2$N, levels = c("0.00 N", "1.10 N", "1.80 N", "2.50 N"))
Tukey_MDA2## N groups
## 1 0.00 N a
## 2 1.10 N b
## 3 2.50 N b
## 4 1.80 N b# Grafico MDA 2
Resumen_MDA2 <- MDA_Agraz2 %>%
dplyr::group_by(., N) %>%
dplyr::summarise( n=n(), media=mean(M2.2), sd=sd(M2.2)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_MDA2, by="N")
head(Resumen_MDA2)## # A tibble: 4 × 7
## N n media sd se ic groups
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 0.00 N 9 67.8 24.0 8.00 18.5 a
## 2 1.10 N 9 15.4 8.86 2.95 6.81 b
## 3 1.80 N 9 13.8 4.74 1.58 3.64 b
## 4 2.50 N 9 13.7 3.99 1.33 3.07 bFigura_MDA2 <- ggplot(Resumen_MDA2, aes(x=N, y=media) )+
geom_errorbar( aes(x=N, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=N, y=media)) +
geom_text( aes(label=groups, x=N, y=media+se+5 ), color="red", size=4)+
labs(x = "Niveles de nitrogeno",
y = (expression("MDA (mmol"~" g"^{-1}~")"))) +
ylim(0.0, 100)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()
Figura_MDA2
### Kruskall MDA
# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
MDA2.TTO <- kruskal.test(MDA_Agraz2$M2.2, MDA_Agraz2$TTO); MDA2.TTO##
## Kruskal-Wallis rank sum test
##
## data: MDA_Agraz2$M2.2 and MDA_Agraz2$TTO
## Kruskal-Wallis chi-squared = 21.572, df = 11, p-value = 0.02791MDA2.N <- kruskal.test(MDA_Agraz2$M2.2, MDA_Agraz2$N); MDA2.N##
## Kruskal-Wallis rank sum test
##
## data: MDA_Agraz2$M2.2 and MDA_Agraz2$N
## Kruskal-Wallis chi-squared = 19.452, df = 3, p-value = 0.0002204MDA2.P <- kruskal.test(MDA_Agraz2$M2.2, MDA_Agraz2$TTO); MDA2.P##
## Kruskal-Wallis rank sum test
##
## data: MDA_Agraz2$M2.2 and MDA_Agraz2$TTO
## Kruskal-Wallis chi-squared = 21.572, df = 11, p-value = 0.02791# ver comparaciones puntuales entre tratamientos
KWlt.MDA2.TTO <- kruskal(MDA_Agraz2$M2.2, MDA_Agraz2$TTO, console=T,
p.adj="bonferroni",group=T)##
## Study: MDA_Agraz2$M2.2 ~ MDA_Agraz2$TTO
## Kruskal-Wallis test's
## Ties or no Ties
##
## Critical Value: 21.57167
## Degrees of freedom: 11
## Pvalue Chisq : 0.0279096
##
## MDA_Agraz2$TTO, means of the ranks
##
## MDA_Agraz2.M2.2 r
## 0.00 N/0.00 P 31.000000 3
## 0.00 N/0.25 P 33.333333 3
## 0.00 N/0.50 P 31.333333 3
## 1.10 N/0.00 P 11.333333 3
## 1.10 N/0.25 P 12.666667 3
## 1.10 N/0.50 P 19.000000 3
## 1.80 N/0.00 P 12.000000 3
## 1.80 N/0.25 P 13.000000 3
## 1.80 N/0.50 P 17.833333 3
## 2.50 N/0.00 P 15.333333 3
## 2.50 N/0.25 P 9.666667 3
## 2.50 N/0.50 P 15.500000 3
##
## Post Hoc Analysis
##
## P value adjustment method: bonferroni
## t-Student: 3.855996
## Alpha : 0.05
## Minimum Significant Difference: 24.78781
##
## Treatments with the same letter are not significantly different.
##
## MDA_Agraz2$M2.2 groups
## 0.00 N/0.25 P 33.333333 a
## 0.00 N/0.50 P 31.333333 a
## 0.00 N/0.00 P 31.000000 a
## 1.10 N/0.50 P 19.000000 a
## 1.80 N/0.50 P 17.833333 a
## 2.50 N/0.50 P 15.500000 a
## 2.50 N/0.00 P 15.333333 a
## 1.80 N/0.25 P 13.000000 a
## 1.10 N/0.25 P 12.666667 a
## 1.80 N/0.00 P 12.000000 a
## 1.10 N/0.00 P 11.333333 a
## 2.50 N/0.25 P 9.666667 aKWlt.MDA2.N <- kruskal(MDA_Agraz2$M2.2, MDA_Agraz2$N, console=T,
p.adj="bonferroni",group=T)##
## Study: MDA_Agraz2$M2.2 ~ MDA_Agraz2$N
## Kruskal-Wallis test's
## Ties or no Ties
##
## Critical Value: 19.45247
## Degrees of freedom: 3
## Pvalue Chisq : 0.00022039
##
## MDA_Agraz2$N, means of the ranks
##
## MDA_Agraz2.M2.2 r
## 0.00 N 31.88889 9
## 1.10 N 14.33333 9
## 1.80 N 14.27778 9
## 2.50 N 13.50000 9
##
## Post Hoc Analysis
##
## P value adjustment method: bonferroni
## t-Student: 2.812343
## Alpha : 0.05
## Minimum Significant Difference: 9.72657
##
## Treatments with the same letter are not significantly different.
##
## MDA_Agraz2$M2.2 groups
## 0.00 N 31.88889 a
## 1.10 N 14.33333 b
## 1.80 N 14.27778 b
## 2.50 N 13.50000 bKWlt.MDA2.P <- kruskal(MDA_Agraz2$M2.2, MDA_Agraz2$P, console=T,
p.adj="bonferroni",group=T)##
## Study: MDA_Agraz2$M2.2 ~ MDA_Agraz2$P
## Kruskal-Wallis test's
## Ties or no Ties
##
## Critical Value: 0.9522888
## Degrees of freedom: 2
## Pvalue Chisq : 0.6211738
##
## MDA_Agraz2$P, means of the ranks
##
## MDA_Agraz2.M2.2 r
## 0.00 P 17.41667 12
## 0.25 P 17.16667 12
## 0.50 P 20.91667 12
##
## Post Hoc Analysis
##
## P value adjustment method: bonferroni
## t-Student: 2.522213
## Alpha : 0.05
## Minimum Significant Difference: 11.00867
##
## Treatments with the same letter are not significantly different.
##
## MDA_Agraz2$M2.2 groups
## 0.50 P 20.91667 a
## 0.00 P 17.41667 a
## 0.25 P 17.16667 aE N S A Y O 2
Cargar base de datos
BD3 <- read_excel("BD3.complete.fisiologicos.xlsx")
BD3.aov <- BD3 %>%
dplyr::select(., -B, -PL, -FC, -R, -CCI, -FC, -CCI.I, -FC.I) %>%
mutate(., TTO=as.factor(TTO), N=as.factor(N), P=as.factor(P), SDT=as.factor(SDT))
str(BD3.aov)## tibble [942 × 7] (S3: tbl_df/tbl/data.frame)
## $ TTO.SDT: chr [1:942] "0.0N/0.0P_6" "0.0N/0.0P_6" "0.0N/0.0P_6" "0.0N/0.0P_6" ...
## $ SDT : Factor w/ 3 levels "16","24","6": 3 3 3 3 3 3 3 3 3 3 ...
## $ TTO : Factor w/ 9 levels "0.0N/0.0P","0.0N/0.25P",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ N : Factor w/ 3 levels "0.0N","1.1N",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ P : Factor w/ 3 levels "0.0P","0.25P",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ CCI.R : num [1:942] 21.4 23.7 18.6 21 22.3 12.6 14.8 18.7 23 26.1 ...
## $ FC.R : num [1:942] 0.794 0.794 0.794 0.794 0.812 0.798 0.826 0.814 0.82 0.826 ...Para CCI
Exploratorio CCI
ggplot(BD3, aes(x=factor(SDT, levels = c("6", "16", "24")), y=CCI.R))+
geom_boxplot(alpha=0.2)+
geom_jitter(aes(color=B), width = 0.1, alpha=0.5)+
facet_grid(~TTO) + ylim(c(0, 55))
24 SDT
# Subconjunto base datos
BD3_24 <- as.data.frame(subset(BD3, subset = SDT==24,
select = c(TTO, N, P, CCI.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_CCI3_24 <- data.frame(lam24 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam24 <- seq(from = -1, to = 20, by = 0.5)
for (i in lam24) {
N <- bartlett.test((CCI.R)^i ~ N, data = BD3_24)
P <- bartlett.test((CCI.R)^i ~ P, data = BD3_24)
# Añadir los valores al dataframe
df_CCI3_24 <- rbind(df_CCI3_24, data.frame(lam24 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_CCI3_24$N.p.value~df_CCI3_24$lam24, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_CCI3_24$P.p.value~df_CCI3_24$lam24, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(CCI.R ~ N * P, data = BD3_24, lambda = seq(-1, 20, 1/10))
#lambda.CCI3_24 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.CCI3_24 <- 0.6
# Creación del modelo ANOVA de dos factores tipo III
mod_CCI.24 <- aov( (CCI.R)^lambda.CCI3_24 ~ N * P, data = BD3_24)
AOV_CCI.24 <- Anova(mod_CCI.24, type = "III"); AOV_CCI.24## Anova Table (Type III tests)
##
## Response: (CCI.R)^lambda.CCI3_24
## Sum Sq Df F value Pr(>F)
## (Intercept) 1136.06 1 358.8405 < 2.2e-16 ***
## N 115.49 2 18.2402 3.398e-08 ***
## P 4.60 2 0.7272 0.4841
## N:P 103.75 4 8.1924 2.815e-06 ***
## Residuals 933.95 295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_CCI.24 <- mod_CCI.24$residuals
res.est_CCI.24 <-res_CCI.24/sd(res_CCI.24)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD3_24$CCI.R), main="Variable original" )
hist( (BD3_24$CCI.R)^lambda.CCI3_24, main="Variable transformada" )
boxplot( (BD3_24$CCI.R), horizontal = T, main="Variable original" )
boxplot( (BD3_24$CCI.R)^lambda.CCI3_24, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_CCI.24, which = 2)
library(nortest)
ad.test(res.est_CCI.24)##
## Anderson-Darling normality test
##
## data: res.est_CCI.24
## A = 3.7685, p-value = 2.043e-09cvm.test(res.est_CCI.24)##
## Cramer-von Mises normality test
##
## data: res.est_CCI.24
## W = 0.64911, p-value = 1.477e-07pearson.test(res.est_CCI.24)##
## Pearson chi-square normality test
##
## data: res.est_CCI.24
## P = 45.605, p-value = 0.0001978ks.test(res.est_CCI.24, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_CCI.24
## D = 0.096275, p-value = 0.007138
## alternative hypothesis: two-sidedshapiro.test(res.est_CCI.24)##
## Shapiro-Wilk normality test
##
## data: res.est_CCI.24
## W = 0.95185, p-value = 1.944e-08# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_CCI.24, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (CCI.R)^lambda.CCI3_24 ~ N, data = BD3_24)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI3_24 by N
## Bartlett's K-squared = 76.024, df = 2, p-value < 2.2e-16bartlett.test( (CCI.R)^lambda.CCI3_24 ~ P, data = BD3_24)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI3_24 by P
## Bartlett's K-squared = 10.763, df = 2, p-value = 0.004602# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_CCI.24 ~ BD3_24$N * BD3_24$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 17.34 < 2.2e-16 ***
## 295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.24 ~ BD3_24$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 23.198 4.269e-10 ***
## 301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.24 ~ BD3_24$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 25.919 4.115e-11 ***
## 301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.24 ~ BD3_24$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 17.34 < 2.2e-16 ***
## 295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test levene - Modelo transformado
leveneTest((BD3_24$CCI.R)^lambda.CCI3_24 ~ BD3_24$N * BD3_24$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 11.313 5.578e-14 ***
## 295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD3_24$CCI.R)^lambda.CCI3_24 ~ BD3_24$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 21.141 2.563e-09 ***
## 301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD3_24$CCI.R)^lambda.CCI3_24 ~ BD3_24$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 7.2008 0.0008816 ***
## 301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.24 ~ BD3_24$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 11.313 5.578e-14 ***
## 295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD3_24$CCI.R)^lambda.CCI3_24 ~ BD3_24$N)##
## hov: Brown-Forsyth
##
## data: (BD3_24$CCI.R)^lambda.CCI3_24
## F = 21.141, df:BD3_24$N = 2, df:Residuals = 301, p-value = 2.563e-09
## alternative hypothesis: variances are not identicalhov((BD3_24$CCI.R)^lambda.CCI3_24 ~ BD3_24$P)##
## hov: Brown-Forsyth
##
## data: (BD3_24$CCI.R)^lambda.CCI3_24
## F = 7.2008, df:BD3_24$P = 2, df:Residuals = 301, p-value = 0.0008816
## alternative hypothesis: variances are not identicalhov((BD3_24$CCI.R)^lambda.CCI3_24 ~ BD3_24$TTO)##
## hov: Brown-Forsyth
##
## data: (BD3_24$CCI.R)^lambda.CCI3_24
## F = 11.313, df:BD3_24$TTO = 8, df:Residuals = 295, p-value = 5.578e-14
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 24 DDT)
library(agricolae)
mod2_CCI.24 <- aov( (CCI.R)^lambda.CCI3_24 ~ TTO, data = BD3_24)
par(mfrow=c(1,2))
plot(mod2_CCI.24, which = 1)
plot(mod2_CCI.24, which = 2)
AOV2_CCI.24 <- Anova(mod2_CCI.24, type = "III")
Tuk_CCI.24 <- HSD.test(mod2_CCI.24, "TTO", console=F, alpha = 0.05)
Tukey_BD3.CCI.24 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_CCI.24$groups[2] ), var = "TTO"),
SDT=rep(x=24, times=9))
Tukey_BD3.CCI.24## TTO groups SDT
## 1 1.8N/0.5P a 24
## 2 1.1N/0.5P ab 24
## 3 1.1N/0.25P b 24
## 4 1.8N/0.25P b 24
## 5 1.8N/0.0P bc 24
## 6 1.1N/0.0P c 24
## 7 0.0N/0.0P d 24
## 8 0.0N/0.25P d 24
## 9 0.0N/0.5P d 24# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD3.CCI.N24 <- kruskal.test(BD3_24$CCI.R, BD3_24$N); KW.BD3.CCI.N24##
## Kruskal-Wallis rank sum test
##
## data: BD3_24$CCI.R and BD3_24$N
## Kruskal-Wallis chi-squared = 124.7, df = 2, p-value < 2.2e-16KW.BD3.CCI.P24 <- kruskal.test(BD3_24$CCI.R, BD3_24$P); KW.BD3.CCI.P24 ##
## Kruskal-Wallis rank sum test
##
## data: BD3_24$CCI.R and BD3_24$P
## Kruskal-Wallis chi-squared = 12.64, df = 2, p-value = 0.0018KW.BD3.CCI.TTO24 <- kruskal.test(BD3_24$CCI.R, BD3_24$TTO); KW.BD3.CCI.TTO24##
## Kruskal-Wallis rank sum test
##
## data: BD3_24$CCI.R and BD3_24$TTO
## Kruskal-Wallis chi-squared = 155.32, df = 8, p-value < 2.2e-16# ver comparaciones puntuales entre tratamientos
KWlt.BD3.CCI.N24 <- kruskal(BD3_24$CCI.R, BD3_24$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.CCI.P24 <- kruskal(BD3_24$CCI.R, BD3_24$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.CCI.TTO24 <- kruskal(BD3_24$CCI.R, BD3_24$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.05)
kruskal_BD3.CCI.24 <- data.frame(
TTO=row.names(KWlt.BD3.CCI.TTO24$groups),
groups_np= c(KWlt.BD3.CCI.TTO24$groups[2])
)
kruskal_BD3.CCI.24## TTO groups
## 1 1.8N/0.5P a
## 2 1.1N/0.5P ab
## 3 1.1N/0.25P b
## 4 1.8N/0.25P bc
## 5 1.8N/0.0P bc
## 6 1.1N/0.0P c
## 7 0.0N/0.25P d
## 8 0.0N/0.0P d
## 9 0.0N/0.5P d# Grafico CCI 24
Resumen_CCI3_24 <- BD3_24 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(CCI.R), media=mean(CCI.R), sd=sd(CCI.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD3.CCI.24, by="TTO") %>%
dplyr::full_join(., kruskal_BD3.CCI.24, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_CCI3_24, n = 12)## # A tibble: 9 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0.0P 36 18 18.0 4.73 0.789 1.60 d 24 d
## 2 0.0N/0.25P 36 18.6 18.1 6.34 1.06 2.15 d 24 d
## 3 0.0N/0.5P 32 15.4 15.4 2.53 0.447 0.913 d 24 d
## 4 1.1N/0.0P 36 28.2 27.9 10.2 1.71 3.47 c 24 c
## 5 1.1N/0.25P 36 36.0 38.0 16.2 2.70 5.49 b 24 b
## 6 1.1N/0.5P 24 45.3 44.2 10.6 2.17 4.50 ab 24 ab
## 7 1.8N/0.0P 32 34.8 34.6 15.1 2.67 5.44 bc 24 bc
## 8 1.8N/0.25P 36 44.2 38.4 18.0 2.99 6.08 b 24 bc
## 9 1.8N/0.5P 36 50.1 50.6 7.16 1.19 2.42 a 24 aFigura_CCI3_24 <- ggplot(Resumen_CCI3_24, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+3.3 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-3.3 ), color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (CCI)")+
ylim(0, 60)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_CCI3_24
### 16 SDT
# Subconjunto base datos
BD3_16 <- as.data.frame(subset(BD3, subset = SDT==16,
select = c(TTO, N, P, CCI.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_CCI3_16 <- data.frame(lam16 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam16 <- seq(from = -10, to = 20, by = 0.5)
for (i in lam16) {
N <- bartlett.test((CCI.R)^i ~ N, data = BD3_16)
P <- bartlett.test((CCI.R)^i ~ P, data = BD3_16)
# Añadir los valores al dataframe
df_CCI3_16 <- rbind(df_CCI3_16, data.frame(lam16 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_CCI3_16$N.p.value~df_CCI3_16$lam16, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_CCI3_16$P.p.value~df_CCI3_16$lam16, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(CCI.R ~ N * P, data = BD3_16, lambda = seq(-10, 20, 1/10))
#lambda.CCI3_16 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.CCI3_16 <- 0.5
# Creación del modelo ANOVA de dos factores tipo III
mod_CCI.16 <- aov( (CCI.R)^lambda.CCI3_16 ~ N * P, data = BD3_16)
AOV_CCI.16 <- Anova(mod_CCI.16, type = "III"); AOV_CCI.16## Anova Table (Type III tests)
##
## Response: (CCI.R)^lambda.CCI3_16
## Sum Sq Df F value Pr(>F)
## (Intercept) 605.63 1 641.0107 < 2.2e-16 ***
## N 29.81 2 15.7777 3.017e-07 ***
## P 3.68 2 1.9485 0.14425
## N:P 9.13 4 2.4148 0.04892 *
## Residuals 288.16 305
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_CCI.16 <- mod_CCI.16$residuals
res.est_CCI.16 <-res_CCI.16/sd(res_CCI.16)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD3_16$CCI.R), main="Variable original" )
hist( (BD3_16$CCI.R)^lambda.CCI3_16, main="Variable transformada" )
boxplot( (BD3_16$CCI.R), horizontal = T, main="Variable original" )
boxplot( (BD3_16$CCI.R)^lambda.CCI3_16, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_CCI.16, which = 2)
library(nortest)
ad.test(res.est_CCI.16)##
## Anderson-Darling normality test
##
## data: res.est_CCI.16
## A = 0.96339, p-value = 0.01496cvm.test(res.est_CCI.16)##
## Cramer-von Mises normality test
##
## data: res.est_CCI.16
## W = 0.1754, p-value = 0.01101pearson.test(res.est_CCI.16)##
## Pearson chi-square normality test
##
## data: res.est_CCI.16
## P = 26, p-value = 0.07446ks.test(res.est_CCI.16, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_CCI.16
## D = 0.057025, p-value = 0.2589
## alternative hypothesis: two-sidedshapiro.test(res.est_CCI.16)##
## Shapiro-Wilk normality test
##
## data: res.est_CCI.16
## W = 0.98782, p-value = 0.009708# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_CCI.16, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (CCI.R)^lambda.CCI3_16 ~ N, data = BD3_16)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI3_16 by N
## Bartlett's K-squared = 9.5414, df = 2, p-value = 0.008474bartlett.test( (CCI.R)^lambda.CCI3_16 ~ P, data = BD3_16)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI3_16 by P
## Bartlett's K-squared = 8.8511, df = 2, p-value = 0.01197# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_CCI.16 ~ BD3_16$N * BD3_16$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 7.4468 4.527e-09 ***
## 305
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.16 ~ BD3_16$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 5.5852 0.004139 **
## 311
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.16 ~ BD3_16$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 2.0406 0.1317
## 311leveneTest(res.est_CCI.16 ~ BD3_16$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 7.4468 4.527e-09 ***
## 305
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test levene - Modelo transformado
leveneTest((BD3_16$CCI.R)^lambda.CCI3_16 ~ BD3_16$N * BD3_16$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 6.8644 2.644e-08 ***
## 305
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD3_16$CCI.R)^lambda.CCI3_16 ~ BD3_16$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 5.6385 0.003932 **
## 311
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD3_16$CCI.R)^lambda.CCI3_16 ~ BD3_16$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 7.2638 0.0008257 ***
## 311
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.16 ~ BD3_16$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 6.8644 2.644e-08 ***
## 305
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD3_16$CCI.R)^lambda.CCI3_16 ~ BD3_16$N)##
## hov: Brown-Forsyth
##
## data: (BD3_16$CCI.R)^lambda.CCI3_16
## F = 5.6385, df:BD3_16$N = 2, df:Residuals = 311, p-value = 0.003932
## alternative hypothesis: variances are not identicalhov((BD3_16$CCI.R)^lambda.CCI3_16 ~ BD3_16$P)##
## hov: Brown-Forsyth
##
## data: (BD3_16$CCI.R)^lambda.CCI3_16
## F = 7.2638, df:BD3_16$P = 2, df:Residuals = 311, p-value = 0.0008257
## alternative hypothesis: variances are not identicalhov((BD3_16$CCI.R)^lambda.CCI3_16 ~ BD3_16$TTO)##
## hov: Brown-Forsyth
##
## data: (BD3_16$CCI.R)^lambda.CCI3_16
## F = 6.8644, df:BD3_16$TTO = 8, df:Residuals = 305, p-value = 2.644e-08
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 16 DDT)
library(agricolae)
mod2_CCI.16 <- aov( (CCI.R)^1 ~ TTO, data = BD3_16)
par(mfrow=c(1,2))
plot(mod2_CCI.16, which = 1)
plot(mod2_CCI.16, which = 2)
AOV2_CCI.16 <- Anova(mod2_CCI.16, type = "III")
Tuk_CCI.16 <- HSD.test(mod2_CCI.16, "TTO", console=F, alpha = 0.05)
Tukey_BD3.CCI.16 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_CCI.16$groups[2] ), var = "TTO"),
SDT=rep(x=16, times=9))
Tukey_BD3.CCI.16## TTO groups SDT
## 1 1.8N/0.5P a 16
## 2 1.1N/0.5P ab 16
## 3 1.8N/0.0P ab 16
## 4 1.1N/0.25P ab 16
## 5 1.8N/0.25P ab 16
## 6 1.1N/0.0P b 16
## 7 0.0N/0.0P c 16
## 8 0.0N/0.25P c 16
## 9 0.0N/0.5P c 16# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD3.CCI.N16 <- kruskal.test(BD3_16$CCI.R, BD3_16$N); KW.BD3.CCI.N16##
## Kruskal-Wallis rank sum test
##
## data: BD3_16$CCI.R and BD3_16$N
## Kruskal-Wallis chi-squared = 115.49, df = 2, p-value < 2.2e-16KW.BD3.CCI.P16 <- kruskal.test(BD3_16$CCI.R, BD3_16$P); KW.BD3.CCI.P16 ##
## Kruskal-Wallis rank sum test
##
## data: BD3_16$CCI.R and BD3_16$P
## Kruskal-Wallis chi-squared = 0.20162, df = 2, p-value = 0.9041KW.BD3.CCI.TTO16 <- kruskal.test(BD3_16$CCI.R, BD3_16$TTO); KW.BD3.CCI.TTO16##
## Kruskal-Wallis rank sum test
##
## data: BD3_16$CCI.R and BD3_16$TTO
## Kruskal-Wallis chi-squared = 120.61, df = 8, p-value < 2.2e-16# ver comparaciones puntuales entre tratamientos
KWlt.BD3.CCI.N16 <- kruskal(BD3_16$CCI.R, BD3_16$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.CCI.P16 <- kruskal(BD3_16$CCI.R, BD3_16$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.CCI.TTO16 <- kruskal(BD3_16$CCI.R, BD3_16$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.05)
kruskal_BD3.CCI.16 <- data.frame(
TTO=row.names(KWlt.BD3.CCI.TTO16$groups),
groups_np= c(KWlt.BD3.CCI.TTO16$groups[2])
)
kruskal_BD3.CCI.16## TTO groups
## 1 1.8N/0.5P a
## 2 1.8N/0.0P a
## 3 1.1N/0.5P a
## 4 1.1N/0.25P a
## 5 1.8N/0.25P a
## 6 1.1N/0.0P a
## 7 0.0N/0.0P b
## 8 0.0N/0.25P b
## 9 0.0N/0.5P b# Grafico CCI 16
Resumen_CCI3_16 <- BD3_16 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(CCI.R), media=mean(CCI.R), sd=sd(CCI.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD3.CCI.16, by="TTO") %>%
dplyr::full_join(., kruskal_BD3.CCI.16, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_CCI3_16, n = 12)## # A tibble: 9 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0.0P 32 21.6 19.5 6.56 1.16 2.36 c 16 b
## 2 0.0N/0.25P 36 16.2 16.8 8.11 1.35 2.75 c 16 b
## 3 0.0N/0.5P 36 15.2 15.8 5.34 0.890 1.81 c 16 b
## 4 1.1N/0.0P 32 26.6 28.4 7.31 1.29 2.63 b 16 a
## 5 1.1N/0.25P 36 32.2 32.9 13.5 2.25 4.58 ab 16 a
## 6 1.1N/0.5P 36 35.4 33.6 14.2 2.36 4.80 ab 16 a
## 7 1.8N/0.0P 35 30.7 33.1 11.7 1.98 4.02 ab 16 a
## 8 1.8N/0.25P 35 31.6 30.2 5.65 0.955 1.94 ab 16 a
## 9 1.8N/0.5P 36 36.0 36.5 12.5 2.09 4.23 a 16 aFigura_CCI3_16 <- ggplot(Resumen_CCI3_16, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+3.3 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-3.3 ), color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (CCI)")+
ylim(0.0, 60.0)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_CCI3_16
06 SDT
# Subconjunto base datos
BD3_06 <- as.data.frame(subset(BD3, subset = SDT==06,
select = c(TTO, N, P, CCI.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_CCI3_06 <- data.frame(lam06 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam06 <- seq(from = -10, to = 20, by = 0.5)
for (i in lam06) {
N <- bartlett.test((CCI.R)^i ~ N, data = BD3_06)
P <- bartlett.test((CCI.R)^i ~ P, data = BD3_06)
# Añadir los valores al dataframe
df_CCI3_06 <- rbind(df_CCI3_06, data.frame(lam06 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_CCI3_06$N.p.value~df_CCI3_06$lam06, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_CCI3_06$P.p.value~df_CCI3_06$lam06, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(CCI.R ~ N * P, data = BD3_06, lambda = seq(-10, 20, 1/10))
#lambda.CCI3_06 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.CCI3_06 <- 0.5
# Creación del modelo ANOVA de dos factores tipo III
mod_CCI.06 <- aov( (CCI.R)^lambda.CCI3_06 ~ N * P, data = BD3_06)
AOV_CCI.06 <- Anova(mod_CCI.06, type = "III"); AOV_CCI.06## Anova Table (Type III tests)
##
## Response: (CCI.R)^lambda.CCI3_06
## Sum Sq Df F value Pr(>F)
## (Intercept) 744.60 1 802.8493 < 2.2e-16 ***
## N 11.97 2 6.4522 0.0017936 **
## P 8.23 2 4.4388 0.0125574 *
## N:P 18.99 4 5.1198 0.0005218 ***
## Residuals 292.14 315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_CCI.06 <- mod_CCI.06$residuals
res.est_CCI.06 <-res_CCI.06/sd(res_CCI.06)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD3_06$CCI.R), main="Variable original" )
hist( (BD3_06$CCI.R)^lambda.CCI3_06, main="Variable transformada" )
boxplot( (BD3_06$CCI.R), horizontal = T, main="Variable original" )
boxplot( (BD3_06$CCI.R)^lambda.CCI3_06, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_CCI.06, which = 2)
library(nortest)
ad.test(res.est_CCI.06)##
## Anderson-Darling normality test
##
## data: res.est_CCI.06
## A = 0.91071, p-value = 0.0202cvm.test(res.est_CCI.06)##
## Cramer-von Mises normality test
##
## data: res.est_CCI.06
## W = 0.16343, p-value = 0.01576pearson.test(res.est_CCI.06)##
## Pearson chi-square normality test
##
## data: res.est_CCI.06
## P = 38.833, p-value = 0.003ks.test(res.est_CCI.06, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_CCI.06
## D = 0.060004, p-value = 0.1938
## alternative hypothesis: two-sidedshapiro.test(res.est_CCI.06)##
## Shapiro-Wilk normality test
##
## data: res.est_CCI.06
## W = 0.9903, p-value = 0.03058# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_CCI.06, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (CCI.R)^lambda.CCI3_06 ~ N, data = BD3_06)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI3_06 by N
## Bartlett's K-squared = 8.7274, df = 2, p-value = 0.01273bartlett.test( (CCI.R)^lambda.CCI3_06 ~ P, data = BD3_06)##
## Bartlett test of homogeneity of variances
##
## data: (CCI.R)^lambda.CCI3_06 by P
## Bartlett's K-squared = 19.145, df = 2, p-value = 6.961e-05# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_CCI.06 ~ BD3_06$N * BD3_06$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 9.559 7.419e-12 ***
## 315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.06 ~ BD3_06$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 4.4592 0.0123 *
## 321
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.06 ~ BD3_06$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 11.318 1.779e-05 ***
## 321
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.06 ~ BD3_06$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 9.559 7.419e-12 ***
## 315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test levene - Modelo transformado
leveneTest((BD3_06$CCI.R)^lambda.CCI3_06 ~ BD3_06$N * BD3_06$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 9.1486 2.522e-11 ***
## 315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD3_06$CCI.R)^lambda.CCI3_06 ~ BD3_06$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 3.2941 0.03836 *
## 321
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD3_06$CCI.R)^lambda.CCI3_06 ~ BD3_06$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 13.787 1.802e-06 ***
## 321
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_CCI.06 ~ BD3_06$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 9.1486 2.522e-11 ***
## 315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD3_06$CCI.R)^lambda.CCI3_06 ~ BD3_06$N)##
## hov: Brown-Forsyth
##
## data: (BD3_06$CCI.R)^lambda.CCI3_06
## F = 3.2941, df:BD3_06$N = 2, df:Residuals = 321, p-value = 0.03836
## alternative hypothesis: variances are not identicalhov((BD3_06$CCI.R)^lambda.CCI3_06 ~ BD3_06$P)##
## hov: Brown-Forsyth
##
## data: (BD3_06$CCI.R)^lambda.CCI3_06
## F = 13.787, df:BD3_06$P = 2, df:Residuals = 321, p-value = 1.802e-06
## alternative hypothesis: variances are not identicalhov((BD3_06$CCI.R)^lambda.CCI3_06 ~ BD3_06$TTO)##
## hov: Brown-Forsyth
##
## data: (BD3_06$CCI.R)^lambda.CCI3_06
## F = 9.1486, df:BD3_06$TTO = 8, df:Residuals = 315, p-value = 2.522e-11
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 06 DDT)
library(agricolae)
mod2_CCI.06 <- aov( (CCI.R)^lambda.CCI3_06 ~ TTO, data = BD3_06)
par(mfrow=c(1,2))
plot(mod2_CCI.06, which = 1)
plot(mod2_CCI.06, which = 2)
AOV2_CCI.06 <- Anova(mod2_CCI.06, type = "III")
Tuk_CCI.06 <- HSD.test(mod2_CCI.06, "TTO", console=F, alpha = 0.05)
Tukey_BD3.CCI.06 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_CCI.06$groups[2] ), var = "TTO"),
SDT=rep(x=06, times=9))
Tukey_BD3.CCI.06## TTO groups SDT
## 1 1.1N/0.0P a 6
## 2 1.1N/0.25P ab 6
## 3 1.8N/0.5P ab 6
## 4 0.0N/0.0P abc 6
## 5 1.1N/0.5P abc 6
## 6 1.8N/0.25P abc 6
## 7 0.0N/0.5P bc 6
## 8 1.8N/0.0P bc 6
## 9 0.0N/0.25P c 6# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD3.CCI.N06 <- kruskal.test(BD3_06$CCI.R, BD3_06$N); KW.BD3.CCI.N06##
## Kruskal-Wallis rank sum test
##
## data: BD3_06$CCI.R and BD3_06$N
## Kruskal-Wallis chi-squared = 13.692, df = 2, p-value = 0.001064KW.BD3.CCI.P06 <- kruskal.test(BD3_06$CCI.R, BD3_06$P); KW.BD3.CCI.P06 ##
## Kruskal-Wallis rank sum test
##
## data: BD3_06$CCI.R and BD3_06$P
## Kruskal-Wallis chi-squared = 1.8399, df = 2, p-value = 0.3985KW.BD3.CCI.TTO06 <- kruskal.test(BD3_06$CCI.R, BD3_06$TTO); KW.BD3.CCI.TTO06##
## Kruskal-Wallis rank sum test
##
## data: BD3_06$CCI.R and BD3_06$TTO
## Kruskal-Wallis chi-squared = 33.552, df = 8, p-value = 4.896e-05# ver comparaciones puntuales entre tratamientos
KWlt.BD3.CCI.N06 <- kruskal(BD3_06$CCI.R, BD3_06$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.CCI.P06 <- kruskal(BD3_06$CCI.R, BD3_06$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.CCI.TTO06 <- kruskal(BD3_06$CCI.R, BD3_06$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.05)
kruskal_BD3.CCI.06 <- data.frame(
TTO=row.names(KWlt.BD3.CCI.TTO06$groups),
groups_np= c(KWlt.BD3.CCI.TTO06$groups[2])
)
kruskal_BD3.CCI.06## TTO groups
## 1 1.8N/0.5P a
## 2 1.1N/0.0P a
## 3 1.1N/0.25P a
## 4 0.0N/0.0P ab
## 5 1.1N/0.5P ab
## 6 1.8N/0.25P ab
## 7 0.0N/0.5P ab
## 8 1.8N/0.0P ab
## 9 0.0N/0.25P b# Grafico CCI 06
Resumen_CCI3_06 <- BD3_06 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(CCI.R), media=mean(CCI.R), sd=sd(CCI.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD3.CCI.06, by="TTO") %>%
dplyr::full_join(., kruskal_BD3.CCI.06, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_CCI3_06, n = 12)## # A tibble: 9 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0.0P 36 21.2 20.9 4.48 0.747 1.52 abc 6 ab
## 2 0.0N/0.25P 36 14.2 15.9 7.75 1.29 2.62 c 6 b
## 3 0.0N/0.5P 36 17.5 18.7 7.17 1.19 2.42 bc 6 ab
## 4 1.1N/0.0P 36 24.4 26.1 12.2 2.04 4.14 a 6 a
## 5 1.1N/0.25P 36 25.6 24.3 9.76 1.63 3.30 ab 6 a
## 6 1.1N/0.5P 36 20.8 20.4 7.05 1.18 2.39 abc 6 ab
## 7 1.8N/0.0P 36 16.8 18.2 8.37 1.39 2.83 bc 6 ab
## 8 1.8N/0.25P 36 18.8 20.0 12.1 2.01 4.09 abc 6 ab
## 9 1.8N/0.5P 36 24.4 23.4 4.94 0.823 1.67 ab 6 aFigura_CCI3_06 <- ggplot(Resumen_CCI3_06, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+3.3 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-3.3 ), color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (CCI)")+
ylim(0.0, 60.0)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_CCI3_06
par(mfrow=c(1,3))
Figura_CCI3_06
Figura_CCI3_16
Figura_CCI3_24
## Para FC ### Exploratorio FC
ggplot(BD3, aes(x=factor(SDT, levels = c("6", "16", "24")), y=FC.R))+
geom_boxplot(alpha=0.2)+
geom_jitter(aes(color=B), width = 0.1, alpha=0.5)+
facet_grid(~TTO) + ylim(c(0, 1))
### 24 SDT
# Subconjunto base datos
BD3_24 <- as.data.frame(subset(BD3, subset = SDT==24,
select = c(TTO, N, P, FC.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_FC3_24 <- data.frame(lam24 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam24 <- seq(from = -1, to = 20, by = 0.5)
for (i in lam24) {
N <- bartlett.test((FC.R)^i ~ N, data = BD3_24)
P <- bartlett.test((FC.R)^i ~ P, data = BD3_24)
# Añadir los valores al dataframe
df_FC3_24 <- rbind(df_FC3_24, data.frame(lam24 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_FC3_24$N.p.value~df_FC3_24$lam24, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_FC3_24$P.p.value~df_FC3_24$lam24, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(FC.R ~ N * P, data = BD3_24, lambda = seq(-1, 20, 1/10))
#lambda.FC3_24 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.FC3_24 <- 5
# Creación del modelo ANOVA de dos factores tipo III
mod_FC.24 <- aov( (FC.R)^lambda.FC3_24 ~ N * P, data = BD3_24)
AOV_FC.24 <- Anova(mod_FC.24, type = "III"); AOV_FC.24## Anova Table (Type III tests)
##
## Response: (FC.R)^lambda.FC3_24
## Sum Sq Df F value Pr(>F)
## (Intercept) 4.1494 1 722.5783 < 2.2e-16 ***
## N 0.0138 2 1.1992 0.302908
## P 0.0214 2 1.8623 0.157136
## N:P 0.0937 4 4.0785 0.003111 **
## Residuals 1.6940 295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_FC.24 <- mod_FC.24$residuals
res.est_FC.24 <-res_FC.24/sd(res_FC.24)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD3_24$FC.R), main="Variable original" )
hist( (BD3_24$FC.R)^lambda.FC3_24, main="Variable transformada" )
boxplot( (BD3_24$FC.R), horizontal = T, main="Variable original" )
boxplot( (BD3_24$FC.R)^lambda.FC3_24, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_FC.24, which = 2)
library(nortest)
ad.test(res.est_FC.24)##
## Anderson-Darling normality test
##
## data: res.est_FC.24
## A = 1.2315, p-value = 0.00326cvm.test(res.est_FC.24)##
## Cramer-von Mises normality test
##
## data: res.est_FC.24
## W = 0.25147, p-value = 0.001231pearson.test(res.est_FC.24)##
## Pearson chi-square normality test
##
## data: res.est_FC.24
## P = 37.579, p-value = 0.002807ks.test(res.est_FC.24, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_FC.24
## D = 0.064316, p-value = 0.1616
## alternative hypothesis: two-sidedshapiro.test(res.est_FC.24)##
## Shapiro-Wilk normality test
##
## data: res.est_FC.24
## W = 0.98886, p-value = 0.01977# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_FC.24, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (FC.R)^lambda.FC3_24 ~ N, data = BD3_24)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC3_24 by N
## Bartlett's K-squared = 20.161, df = 2, p-value = 4.188e-05bartlett.test( (FC.R)^lambda.FC3_24 ~ P, data = BD3_24)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC3_24 by P
## Bartlett's K-squared = 0.42109, df = 2, p-value = 0.8101# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_FC.24 ~ BD3_24$N * BD3_24$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 4.0269 0.0001475 ***
## 295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.24 ~ BD3_24$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 8.7151 0.0002092 ***
## 301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.24 ~ BD3_24$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 0.2769 0.7583
## 301leveneTest(res.est_FC.24 ~ BD3_24$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 4.0269 0.0001475 ***
## 295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test levene - Modelo transformado
leveneTest((BD3_24$FC.R)^lambda.FC3_24 ~ BD3_24$N * BD3_24$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 3.3187 0.001192 **
## 295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD3_24$FC.R)^lambda.FC3_24 ~ BD3_24$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 10.077 5.805e-05 ***
## 301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD3_24$FC.R)^lambda.FC3_24 ~ BD3_24$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 0.3148 0.7302
## 301leveneTest(res.est_FC.24 ~ BD3_24$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 3.3187 0.001192 **
## 295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD3_24$FC.R)^lambda.FC3_24 ~ BD3_24$N)##
## hov: Brown-Forsyth
##
## data: (BD3_24$FC.R)^lambda.FC3_24
## F = 10.077, df:BD3_24$N = 2, df:Residuals = 301, p-value = 5.805e-05
## alternative hypothesis: variances are not identicalhov((BD3_24$FC.R)^lambda.FC3_24 ~ BD3_24$P)##
## hov: Brown-Forsyth
##
## data: (BD3_24$FC.R)^lambda.FC3_24
## F = 0.31478, df:BD3_24$P = 2, df:Residuals = 301, p-value = 0.7302
## alternative hypothesis: variances are not identicalhov((BD3_24$FC.R)^lambda.FC3_24 ~ BD3_24$TTO)##
## hov: Brown-Forsyth
##
## data: (BD3_24$FC.R)^lambda.FC3_24
## F = 3.3187, df:BD3_24$TTO = 8, df:Residuals = 295, p-value = 0.001192
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 24 DDT)
library(agricolae)
mod2_FC.24 <- aov( (FC.R)^lambda.FC3_24 ~ TTO, data = BD3_24)
par(mfrow=c(1,2))
plot(mod2_FC.24, which = 1)
plot(mod2_FC.24, which = 2)
AOV2_FC.24 <- Anova(mod2_FC.24, type = "III")
Tuk_FC.24 <- HSD.test(mod2_FC.24, "TTO", console=F, alpha = 0.05)
Tukey_BD3.FC.24 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_FC.24$groups[2] ), var = "TTO"),
SDT=rep(x=24, times=9))
Tukey_BD3.FC.24## TTO groups SDT
## 1 1.1N/0.25P a 24
## 2 1.1N/0.5P ab 24
## 3 1.8N/0.25P abc 24
## 4 1.8N/0.5P abc 24
## 5 0.0N/0.5P abc 24
## 6 0.0N/0.0P bc 24
## 7 0.0N/0.25P c 24
## 8 1.1N/0.0P c 24
## 9 1.8N/0.0P c 24# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD3.FC.N24 <- kruskal.test(BD3_24$FC.R, BD3_24$N); KW.BD3.FC.N24##
## Kruskal-Wallis rank sum test
##
## data: BD3_24$FC.R and BD3_24$N
## Kruskal-Wallis chi-squared = 8.556, df = 2, p-value = 0.01387KW.BD3.FC.P24 <- kruskal.test(BD3_24$FC.R, BD3_24$P); KW.BD3.FC.P24 ##
## Kruskal-Wallis rank sum test
##
## data: BD3_24$FC.R and BD3_24$P
## Kruskal-Wallis chi-squared = 22.815, df = 2, p-value = 1.111e-05KW.BD3.FC.TTO24 <- kruskal.test(BD3_24$FC.R, BD3_24$TTO); KW.BD3.FC.TTO24##
## Kruskal-Wallis rank sum test
##
## data: BD3_24$FC.R and BD3_24$TTO
## Kruskal-Wallis chi-squared = 47.496, df = 8, p-value = 1.233e-07# ver comparaciones puntuales entre tratamientos
KWlt.BD3.FC.N24 <- kruskal(BD3_24$FC.R, BD3_24$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.FC.P24 <- kruskal(BD3_24$FC.R, BD3_24$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.FC.TTO24 <- kruskal(BD3_24$FC.R, BD3_24$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.05)
kruskal_BD3.FC.24 <- data.frame(
TTO=row.names(KWlt.BD3.FC.TTO24$groups),
groups_np= c(KWlt.BD3.FC.TTO24$groups[2])
)
kruskal_BD3.FC.24## TTO groups
## 1 1.1N/0.25P a
## 2 1.1N/0.5P ab
## 3 1.8N/0.25P abc
## 4 1.8N/0.5P abc
## 5 0.0N/0.5P abc
## 6 0.0N/0.0P bc
## 7 0.0N/0.25P c
## 8 1.1N/0.0P c
## 9 1.8N/0.0P c# Grafico FC 24
Resumen_FC3_24 <- BD3_24 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(FC.R), media=mean(FC.R), sd=sd(FC.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD3.FC.24, by="TTO") %>%
dplyr::full_join(., kruskal_BD3.FC.24, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_FC3_24, n = 12)## # A tibble: 9 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0.… 36 0.804 0.805 0.0192 0.00320 0.00649 bc 24 bc
## 2 0.0N/0.… 36 0.803 0.793 0.0391 0.00651 0.0132 c 24 c
## 3 0.0N/0.… 32 0.810 0.812 0.0297 0.00526 0.0107 abc 24 abc
## 4 1.1N/0.… 36 0.786 0.792 0.0388 0.00646 0.0131 c 24 c
## 5 1.1N/0.… 36 0.832 0.831 0.0214 0.00357 0.00725 a 24 a
## 6 1.1N/0.… 24 0.834 0.829 0.0295 0.00603 0.0125 ab 24 ab
## 7 1.8N/0.… 32 0.781 0.786 0.0500 0.00884 0.0180 c 24 c
## 8 1.8N/0.… 36 0.818 0.811 0.0471 0.00786 0.0159 abc 24 abc
## 9 1.8N/0.… 36 0.818 0.810 0.0433 0.00721 0.0146 abc 24 abcFigura_FC3_24 <- ggplot(Resumen_FC3_24, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+0.06 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-0.06 ), color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (FC)")+
ylim(0.2, 1.0)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_FC3_24
### 16 SDT
# Subconjunto base datos
BD3_16 <- as.data.frame(subset(BD3, subset = SDT==16,
select = c(TTO, N, P, FC.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_FC3_16 <- data.frame(lam16 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam16 <- seq(from = -10, to = 20, by = 0.5)
for (i in lam16) {
N <- bartlett.test((FC.R)^i ~ N, data = BD3_16)
P <- bartlett.test((FC.R)^i ~ P, data = BD3_16)
# Añadir los valores al dataframe
df_FC3_16 <- rbind(df_FC3_16, data.frame(lam16 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_FC3_16$N.p.value~df_FC3_16$lam16, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_FC3_16$P.p.value~df_FC3_16$lam16, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(FC.R ~ N * P, data = BD3_16, lambda = seq(-10, 20, 1/10))
#lambda.FC3_16 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.FC3_16 <- 6.3
# Creación del modelo ANOVA de dos factores tipo III
mod_FC.16 <- aov( (FC.R)^lambda.FC3_16 ~ N * P, data = BD3_16)
AOV_FC.16 <- Anova(mod_FC.16, type = "III"); AOV_FC.16## Anova Table (Type III tests)
##
## Response: (FC.R)^lambda.FC3_16
## Sum Sq Df F value Pr(>F)
## (Intercept) 1.44787 1 518.5000 < 2.2e-16 ***
## N 0.01735 2 3.1061 0.0462 *
## P 0.00407 2 0.7286 0.4834
## N:P 0.07506 4 6.7196 3.403e-05 ***
## Residuals 0.85169 305
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_FC.16 <- mod_FC.16$residuals
res.est_FC.16 <-res_FC.16/sd(res_FC.16)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD3_16$FC.R), main="Variable original" )
hist( (BD3_16$FC.R)^lambda.FC3_16, main="Variable transformada" )
boxplot( (BD3_16$FC.R), horizontal = T, main="Variable original" )
boxplot( (BD3_16$FC.R)^lambda.FC3_16, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_FC.16, which = 2)
library(nortest)
ad.test(res.est_FC.16)##
## Anderson-Darling normality test
##
## data: res.est_FC.16
## A = 1.287, p-value = 0.00238cvm.test(res.est_FC.16)##
## Cramer-von Mises normality test
##
## data: res.est_FC.16
## W = 0.23133, p-value = 0.002168pearson.test(res.est_FC.16)##
## Pearson chi-square normality test
##
## data: res.est_FC.16
## P = 42.688, p-value = 0.0005331ks.test(res.est_FC.16, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_FC.16
## D = 0.06204, p-value = 0.1782
## alternative hypothesis: two-sidedshapiro.test(res.est_FC.16)##
## Shapiro-Wilk normality test
##
## data: res.est_FC.16
## W = 0.98847, p-value = 0.0136# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_FC.16, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (FC.R)^lambda.FC3_16 ~ N, data = BD3_16)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC3_16 by N
## Bartlett's K-squared = 3.057, df = 2, p-value = 0.2169bartlett.test( (FC.R)^lambda.FC3_16 ~ P, data = BD3_16)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC3_16 by P
## Bartlett's K-squared = 13.293, df = 2, p-value = 0.001299# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_FC.16 ~ BD3_16$N * BD3_16$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 1.2867 0.2497
## 305leveneTest(res.est_FC.16 ~ BD3_16$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 0.0105 0.9895
## 311leveneTest(res.est_FC.16 ~ BD3_16$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 2.8726 0.05805 .
## 311
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.16 ~ BD3_16$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 1.2867 0.2497
## 305# Test levene - Modelo transformado
leveneTest((BD3_16$FC.R)^lambda.FC3_16 ~ BD3_16$N * BD3_16$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 1.3696 0.2092
## 305leveneTest((BD3_16$FC.R)^lambda.FC3_16 ~ BD3_16$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 2.2221 0.1101
## 311leveneTest((BD3_16$FC.R)^lambda.FC3_16 ~ BD3_16$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 7.5833 0.0006087 ***
## 311
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.16 ~ BD3_16$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 1.3696 0.2092
## 305# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD3_16$FC.R)^lambda.FC3_16 ~ BD3_16$N)##
## hov: Brown-Forsyth
##
## data: (BD3_16$FC.R)^lambda.FC3_16
## F = 2.2221, df:BD3_16$N = 2, df:Residuals = 311, p-value = 0.1101
## alternative hypothesis: variances are not identicalhov((BD3_16$FC.R)^lambda.FC3_16 ~ BD3_16$P)##
## hov: Brown-Forsyth
##
## data: (BD3_16$FC.R)^lambda.FC3_16
## F = 7.5833, df:BD3_16$P = 2, df:Residuals = 311, p-value = 0.0006087
## alternative hypothesis: variances are not identicalhov((BD3_16$FC.R)^lambda.FC3_16 ~ BD3_16$TTO)##
## hov: Brown-Forsyth
##
## data: (BD3_16$FC.R)^lambda.FC3_16
## F = 1.3696, df:BD3_16$TTO = 8, df:Residuals = 305, p-value = 0.2092
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 16 DDT)
library(agricolae)
mod2_FC.16 <- aov( (FC.R)^lambda.FC3_16 ~ TTO, data = BD3_16)
par(mfrow=c(1,2))
plot(mod2_FC.16, which = 1)
plot(mod2_FC.16, which = 2)
AOV2_FC.16 <- Anova(mod2_FC.16, type = "III")
Tuk_FC.16 <- HSD.test(mod2_FC.16, "TTO", console=F, alpha = 0.05)
Tukey_BD3.FC.16 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_FC.16$groups[2] ), var = "TTO"),
SDT=rep(x=16, times=9))
Tukey_BD3.FC.16## TTO groups SDT
## 1 1.1N/0.25P a 16
## 2 1.1N/0.5P ab 16
## 3 1.8N/0.25P bc 16
## 4 1.8N/0.0P bcd 16
## 5 1.8N/0.5P bcd 16
## 6 1.1N/0.0P bcde 16
## 7 0.0N/0.0P cde 16
## 8 0.0N/0.5P de 16
## 9 0.0N/0.25P e 16# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD3.FC.N16 <- kruskal.test(BD3_16$FC.R, BD3_16$N); KW.BD3.FC.N16##
## Kruskal-Wallis rank sum test
##
## data: BD3_16$FC.R and BD3_16$N
## Kruskal-Wallis chi-squared = 46.882, df = 2, p-value = 6.602e-11KW.BD3.FC.P16 <- kruskal.test(BD3_16$FC.R, BD3_16$P); KW.BD3.FC.P16 ##
## Kruskal-Wallis rank sum test
##
## data: BD3_16$FC.R and BD3_16$P
## Kruskal-Wallis chi-squared = 8.0373, df = 2, p-value = 0.01798KW.BD3.FC.TTO16 <- kruskal.test(BD3_16$FC.R, BD3_16$TTO); KW.BD3.FC.TTO16##
## Kruskal-Wallis rank sum test
##
## data: BD3_16$FC.R and BD3_16$TTO
## Kruskal-Wallis chi-squared = 74.38, df = 8, p-value = 6.564e-13# ver comparaciones puntuales entre tratamientos
KWlt.BD3.FC.N16 <- kruskal(BD3_16$FC.R, BD3_16$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.FC.P16 <- kruskal(BD3_16$FC.R, BD3_16$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.FC.TTO16 <- kruskal(BD3_16$FC.R, BD3_16$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.05)
kruskal_BD3.FC.16 <- data.frame(
TTO=row.names(KWlt.BD3.FC.TTO16$groups),
groups_np= c(KWlt.BD3.FC.TTO16$groups[2])
)
kruskal_BD3.FC.16## TTO groups
## 1 1.1N/0.25P a
## 2 1.1N/0.5P ab
## 3 1.8N/0.25P bc
## 4 1.8N/0.0P bcd
## 5 1.8N/0.5P bcde
## 6 1.1N/0.0P cde
## 7 0.0N/0.0P de
## 8 0.0N/0.5P e
## 9 0.0N/0.25P e# Grafico FC 16
Resumen_FC3_16 <- BD3_16 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(FC.R), media=mean(FC.R), sd=sd(FC.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD3.FC.16, by="TTO") %>%
dplyr::full_join(., kruskal_BD3.FC.16, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_FC3_16, n = 12)## # A tibble: 9 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0.… 32 0.782 0.778 0.0350 0.00619 0.0126 cde 16 de
## 2 0.0N/0.… 36 0.774 0.768 0.0425 0.00709 0.0144 e 16 e
## 3 0.0N/0.… 36 0.783 0.780 0.0245 0.00409 0.00831 de 16 e
## 4 1.1N/0.… 32 0.784 0.785 0.0305 0.00539 0.0110 bcde 16 cde
## 5 1.1N/0.… 36 0.824 0.824 0.0233 0.00389 0.00790 a 16 a
## 6 1.1N/0.… 36 0.808 0.806 0.0267 0.00445 0.00903 ab 16 ab
## 7 1.8N/0.… 35 0.804 0.797 0.0285 0.00481 0.00978 bcd 16 bcd
## 8 1.8N/0.… 35 0.801 0.801 0.0305 0.00516 0.0105 bc 16 bc
## 9 1.8N/0.… 36 0.792 0.798 0.0238 0.00397 0.00806 bcd 16 bcdeFigura_FC3_16 <- ggplot(Resumen_FC3_16, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+0.06 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-0.06 ), color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (FC)")+
ylim(0.2, 1.0)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_FC3_16
### 06 SDT
# Subconjunto base datos
BD3_06 <- as.data.frame(subset(BD3, subset = SDT==06,
select = c(TTO, N, P, FC.R) ) )# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_FC3_06 <- data.frame(lam06 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam06 <- seq(from = -10, to = 20, by = 0.5)
for (i in lam06) {
N <- bartlett.test((FC.R)^i ~ N, data = BD3_06)
P <- bartlett.test((FC.R)^i ~ P, data = BD3_06)
# Añadir los valores al dataframe
df_FC3_06 <- rbind(df_FC3_06, data.frame(lam06 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_FC3_06$N.p.value~df_FC3_06$lam06, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_FC3_06$P.p.value~df_FC3_06$lam06, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
#b <- boxcox(FC.R ~ N * P, data = BD3_06, lambda = seq(-10, 20, 1/10))
#lambda.FC3_06 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.FC3_06 <- 12
# Creación del modelo ANOVA de dos factores tipo III
mod_FC.06 <- aov( (FC.R)^lambda.FC3_06 ~ N * P, data = BD3_06)
AOV_FC.06 <- Anova(mod_FC.06, type = "III"); AOV_FC.06## Anova Table (Type III tests)
##
## Response: (FC.R)^lambda.FC3_06
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.168180 1 262.1684 < 2.2e-16 ***
## N 0.003443 2 2.6832 0.0699032 .
## P 0.000911 2 0.7097 0.4925532
## N:P 0.015287 4 5.9576 0.0001241 ***
## Residuals 0.202071 315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_FC.06 <- mod_FC.06$residuals
res.est_FC.06 <-res_FC.06/sd(res_FC.06)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (BD3_06$FC.R), main="Variable original" )
hist( (BD3_06$FC.R)^lambda.FC3_06, main="Variable transformada" )
boxplot( (BD3_06$FC.R), horizontal = T, main="Variable original" )
boxplot( (BD3_06$FC.R)^lambda.FC3_06, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_FC.06, which = 2)
library(nortest)
ad.test(res.est_FC.06)##
## Anderson-Darling normality test
##
## data: res.est_FC.06
## A = 1.1691, p-value = 0.004651cvm.test(res.est_FC.06)##
## Cramer-von Mises normality test
##
## data: res.est_FC.06
## W = 0.2237, p-value = 0.002694pearson.test(res.est_FC.06)##
## Pearson chi-square normality test
##
## data: res.est_FC.06
## P = 71.111, p-value = 2.929e-08ks.test(res.est_FC.06, "pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: res.est_FC.06
## D = 0.075897, p-value = 0.04785
## alternative hypothesis: two-sidedshapiro.test(res.est_FC.06)##
## Shapiro-Wilk normality test
##
## data: res.est_FC.06
## W = 0.97276, p-value = 8.391e-06# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_FC.06, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (FC.R)^lambda.FC3_06 ~ N, data = BD3_06)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC3_06 by N
## Bartlett's K-squared = 8.6554, df = 2, p-value = 0.0132bartlett.test( (FC.R)^lambda.FC3_06 ~ P, data = BD3_06)##
## Bartlett test of homogeneity of variances
##
## data: (FC.R)^lambda.FC3_06 by P
## Bartlett's K-squared = 7.5397, df = 2, p-value = 0.02306# Tests sin certeza de normalidad (Levene, BrownForsyth)
library(car)
# Test levene - Modelo sin transformar
leveneTest(res.est_FC.06 ~ BD3_06$N * BD3_06$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 4.0672 0.0001265 ***
## 315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.06 ~ BD3_06$N, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 0.6261 0.5353
## 321leveneTest(res.est_FC.06 ~ BD3_06$P, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 3.3764 0.03539 *
## 321
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(res.est_FC.06 ~ BD3_06$TTO, center="mean")## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 8 4.0672 0.0001265 ***
## 315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test levene - Modelo transformado
leveneTest((BD3_06$FC.R)^lambda.FC3_06 ~ BD3_06$N * BD3_06$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 3.7006 0.0003784 ***
## 315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest((BD3_06$FC.R)^lambda.FC3_06 ~ BD3_06$N, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 0.9406 0.3914
## 321leveneTest((BD3_06$FC.R)^lambda.FC3_06 ~ BD3_06$P, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 2 2.0616 0.1289
## 321leveneTest(res.est_FC.06 ~ BD3_06$TTO, center="median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 8 3.7006 0.0003784 ***
## 315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test BrownForsyth - Modelo transformado
library(HH)
hov((BD3_06$FC.R)^lambda.FC3_06 ~ BD3_06$N)##
## hov: Brown-Forsyth
##
## data: (BD3_06$FC.R)^lambda.FC3_06
## F = 0.94065, df:BD3_06$N = 2, df:Residuals = 321, p-value = 0.3914
## alternative hypothesis: variances are not identicalhov((BD3_06$FC.R)^lambda.FC3_06 ~ BD3_06$P)##
## hov: Brown-Forsyth
##
## data: (BD3_06$FC.R)^lambda.FC3_06
## F = 2.0616, df:BD3_06$P = 2, df:Residuals = 321, p-value = 0.1289
## alternative hypothesis: variances are not identicalhov((BD3_06$FC.R)^lambda.FC3_06 ~ BD3_06$TTO)##
## hov: Brown-Forsyth
##
## data: (BD3_06$FC.R)^lambda.FC3_06
## F = 3.7006, df:BD3_06$TTO = 8, df:Residuals = 315, p-value = 0.0003784
## alternative hypothesis: variances are not identical# Analisis parametrico (Tukey 06 DDT)
library(agricolae)
mod2_FC.06 <- aov( (FC.R)^lambda.FC3_06 ~ TTO, data = BD3_06)
par(mfrow=c(1,2))
plot(mod2_FC.06, which = 1)
plot(mod2_FC.06, which = 2)
AOV2_FC.06 <- Anova(mod2_FC.06, type = "III")
Tuk_FC.06 <- HSD.test(mod2_FC.06, "TTO", console=F, alpha = 0.05)
Tukey_BD3.FC.06 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_FC.06$groups[2] ), var = "TTO"),
SDT=rep(x=06, times=9))
Tukey_BD3.FC.06## TTO groups SDT
## 1 1.8N/0.25P a 6
## 2 1.1N/0.25P ab 6
## 3 1.1N/0.5P abc 6
## 4 1.8N/0.0P abc 6
## 5 1.1N/0.0P bcd 6
## 6 0.0N/0.25P bcd 6
## 7 0.0N/0.5P bcd 6
## 8 0.0N/0.0P cd 6
## 9 1.8N/0.5P d 6# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
KW.BD3.FC.N06 <- kruskal.test(BD3_06$FC.R, BD3_06$N); KW.BD3.FC.N06##
## Kruskal-Wallis rank sum test
##
## data: BD3_06$FC.R and BD3_06$N
## Kruskal-Wallis chi-squared = 8.8452, df = 2, p-value = 0.012KW.BD3.FC.P06 <- kruskal.test(BD3_06$FC.R, BD3_06$P); KW.BD3.FC.P06 ##
## Kruskal-Wallis rank sum test
##
## data: BD3_06$FC.R and BD3_06$P
## Kruskal-Wallis chi-squared = 12.059, df = 2, p-value = 0.002407KW.BD3.FC.TTO06 <- kruskal.test(BD3_06$FC.R, BD3_06$TTO); KW.BD3.FC.TTO06##
## Kruskal-Wallis rank sum test
##
## data: BD3_06$FC.R and BD3_06$TTO
## Kruskal-Wallis chi-squared = 41.009, df = 8, p-value = 2.077e-06# ver comparaciones puntuales entre tratamientos
KWlt.BD3.FC.N06 <- kruskal(BD3_06$FC.R, BD3_06$N, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.FC.P06 <- kruskal(BD3_06$FC.R, BD3_06$P, console=F,
p.adj="bonferroni",group=T)
KWlt.BD3.FC.TTO06 <- kruskal(BD3_06$FC.R, BD3_06$TTO, console=F,
p.adj="bonferroni",group=T, alpha = 0.05)
kruskal_BD3.FC.06 <- data.frame(
TTO=row.names(KWlt.BD3.FC.TTO06$groups),
groups_np= c(KWlt.BD3.FC.TTO06$groups[2])
)
kruskal_BD3.FC.06## TTO groups
## 1 1.8N/0.25P a
## 2 1.1N/0.25P ab
## 3 1.1N/0.5P ab
## 4 1.8N/0.0P ab
## 5 1.1N/0.0P abc
## 6 0.0N/0.25P abc
## 7 0.0N/0.5P abc
## 8 0.0N/0.0P bc
## 9 1.8N/0.5P c# Grafico FC 06
Resumen_FC3_06 <- BD3_06 %>%
dplyr::group_by(., TTO) %>%
dplyr::summarise( n=n(), median=median(FC.R), media=mean(FC.R), sd=sd(FC.R)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_BD3.FC.06, by="TTO") %>%
dplyr::full_join(., kruskal_BD3.FC.06, by="TTO", suffix = c(".tuk", ".kru"))
head(Resumen_FC3_06, n = 12)## # A tibble: 9 × 10
## TTO n median media sd se ic groups.tuk SDT groups.kru
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr>
## 1 0.0N/0.… 36 0.794 0.795 0.0290 0.00484 0.00982 cd 6 bc
## 2 0.0N/0.… 36 0.804 0.803 0.0198 0.00330 0.00669 bcd 6 abc
## 3 0.0N/0.… 36 0.807 0.802 0.0202 0.00336 0.00682 bcd 6 abc
## 4 1.1N/0.… 36 0.809 0.799 0.0372 0.00620 0.0126 bcd 6 abc
## 5 1.1N/0.… 36 0.814 0.815 0.0126 0.00210 0.00427 ab 6 ab
## 6 1.1N/0.… 36 0.814 0.813 0.0179 0.00299 0.00607 abc 6 ab
## 7 1.8N/0.… 36 0.807 0.809 0.0204 0.00340 0.00689 abc 6 ab
## 8 1.8N/0.… 36 0.819 0.820 0.0212 0.00353 0.00716 a 6 a
## 9 1.8N/0.… 36 0.789 0.790 0.0237 0.00395 0.00801 d 6 cFigura_FC3_06 <- ggplot(Resumen_FC3_06, aes(x=TTO, y=media) )+
geom_errorbar( aes(x=TTO, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=media)) +
geom_text( aes(label=groups.kru, x=TTO, y=media+se+0.06 ), color="red", size=4)+
#geom_text( aes(label=groups.tuk, x=TTO, y=media-se-0.06 ), color="blue", size=4)+
labs(x = "Tratamientos",
y = "Contenido relativo de clorofila (FC)")+
ylim(0.2, 1.0)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()+
coord_flip()
Figura_FC3_06
par(mfrow=c(1,3))
Figura_FC3_06
Figura_FC3_16
Figura_FC3_24
# Resumen
CCI2 <- c(KW.BD2.CCI.N25$p.value, KW.BD2.CCI.N16$p.value, KW.BD2.CCI.N07$p.value,
KW.BD2.CCI.P25$p.value, KW.BD2.CCI.P16$p.value, KW.BD2.CCI.P07$p.value,
KW.BD2.CCI.TTO25$p.value, KW.BD2.CCI.TTO16$p.value, KW.BD2.CCI.TTO07$p.value)
CCI3 <- c(KW.BD3.CCI.N24$p.value, KW.BD3.CCI.N16$p.value, KW.BD3.CCI.N06$p.value,
KW.BD3.CCI.P24$p.value, KW.BD3.CCI.P16$p.value, KW.BD3.CCI.P06$p.value,
KW.BD3.CCI.TTO24$p.value, KW.BD3.CCI.TTO16$p.value, KW.BD3.CCI.TTO06$p.value)
resumen.KW.CCI <- round(matrix(c(CCI2, CCI3), nrow = 3, ncol = 6, byrow = TRUE), digits = 4)
rownames(resumen.KW.CCI) <- c("N","P","NxP")
colnames(resumen.KW.CCI) <- c("E2-25","E2-16","E2-07"
,"E3-24","E3-25","E3-06")
resumen.KW.CCI## E2-25 E2-16 E2-07 E3-24 E3-25 E3-06
## N 0.0000 0.0000 0.0043 0.1002 0.3201 0.0095
## P 0.0000 0.0000 0.0005 0.0000 0.0000 0.0011
## NxP 0.0018 0.9041 0.3985 0.0000 0.0000 0.0000FC2 <- c(KW.BD2.FC.N25$p.value, KW.BD2.FC.N16$p.value, KW.BD2.FC.N07$p.value,
KW.BD2.FC.P25$p.value, KW.BD2.FC.P16$p.value, KW.BD2.FC.P07$p.value,
KW.BD2.FC.TTO25$p.value, KW.BD2.FC.TTO16$p.value, KW.BD2.FC.TTO07$p.value)
FC3 <- c(KW.BD3.FC.N24$p.value, KW.BD3.FC.N16$p.value, KW.BD3.FC.N06$p.value,
KW.BD3.FC.P24$p.value, KW.BD3.FC.P16$p.value, KW.BD3.FC.P06$p.value,
KW.BD3.FC.TTO24$p.value, KW.BD3.FC.TTO16$p.value, KW.BD3.FC.TTO06$p.value)
resumen.KW.FC <- round(matrix(c(FC2, FC3), nrow = 3, ncol = 6, byrow = TRUE), digits = 4)
rownames(resumen.KW.FC) <- c("N","P","NxP")
colnames(resumen.KW.FC) <- c("E2-25","E2-16","E2-07"
,"E3-24","E3-25","E3-06")
resumen.KW.FC## E2-25 E2-16 E2-07 E3-24 E3-25 E3-06
## N 0 0.0027 0.0001 0.5246 0.0763 0.7265
## P 0 0.0121 0.0000 0.0139 0.0000 0.0120
## NxP 0 0.0180 0.0024 0.0000 0.0000 0.0000M D A 3
Exploratorio MDA
MDA_Agraz3 <- read_excel("MDA Agraz.xlsx", sheet = "MDA ENSAYO 2")
head(MDA_Agraz3)## # A tibble: 6 × 22
## SDT TTO N P PL M V `440-` `532-` `600-` `440+` `532+`
## <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 24 1 0.00 N 0.00… 6 0.242 5 -0.02 -0.036 -0.093 0.206 0.164
## 2 24 1 0.00 N 0.00… 108 0.177 5 0.016 -0.009 -0.02 0.322 0.246
## 3 24 1 0.00 N 0.00… 201 0.179 5 0.149 0.11 0.006 0.706 0.809
## 4 24 2 1.10 N 0.00… 220 0.208 5 0.115 0.05 0.026 0.375 0.286
## 5 24 2 1.10 N 0.00… 114 0.273 5 0.008 -0.013 -0.035 0.373 0.284
## 6 24 2 1.10 N 0.00… NA NA NA NA NA NA NA NA
## # ℹ 10 more variables: `600+` <dbl>, A <dbl>, B <dbl>, C <dbl>, M1 <dbl>,
## # M1.1 <dbl>, M2 <dbl>, M2.1 <dbl>, M2.2 <dbl>, x <dbl>MDA_Agraz3 <- MDA_Agraz3 %>%
unite(., col="TTO", c("N","P"), sep = "/", remove = FALSE)
Resumen3 <- MDA_Agraz3 %>%
dplyr::group_by(., N) %>%
dplyr::summarise(.,
mediaM2=mean(M2.2), medianaM2=median(M2.2), stdM2=sd(M2.2), n2= n(),
mediaM1=mean(M1.1), medianaM1=median(M1.1), stdM1=sd(M1.1), n1= n() )
Resumen3## # A tibble: 3 × 9
## N mediaM2 medianaM2 stdM2 n2 mediaM1 medianaM1 stdM1 n1
## <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <int>
## 1 0.00 N 62.1 53.2 32.5 9 43.0 42.8 34.2 9
## 2 1.10 N 17.0 16 10.1 9 13.0 14 8.61 9
## 3 1.80 N 19.0 18.5 6.23 9 12.7 12 5.23 9ggplot(MDA_Agraz3, aes(x=factor(TTO), y=M2.2) )+
geom_point()+
geom_col(alpha=0.2)+
coord_flip()
ggplot(MDA_Agraz3, aes(x=factor(N), y=M2.2) )+
geom_point()+
geom_col(alpha=0.2)+
coord_flip()
ggplot(MDA_Agraz3, aes(x=factor(P), y=M2.2) )+
geom_point()+
geom_col(alpha=0.2)+
coord_flip()
# Transformacion de la variable para cumplir supuestos
# Mediante iteraciones
df_MDA3 <- data.frame(lam06 = numeric(), N.p.value = numeric(), P.p.value = numeric())
lam3 <- seq(from = -10, to = 20, by = 0.5)
for (i in lam3) {
N <- bartlett.test((M2.2)^i ~ N, data = MDA_Agraz3)
P <- bartlett.test((M2.2)^i ~ P, data = MDA_Agraz3)
# Añadir los valores al dataframe
df_MDA3 <- rbind(df_MDA3, data.frame(lam3 = i,
N.p.value = N$p.value,
P.p.value = P$p.value))
}
par(mfrow=c(1,3))
plot(df_MDA3$N.p.value~df_MDA3$lam3, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
plot(df_MDA3$P.p.value~df_MDA3$lam3, ylim = c(0, 0.5))
abline(h = 0.05, col="red")
# Mediante BOXCOX
library(car)
library(EnvStats)
# b <- boxcox(M2.2 ~ N * P, data = MDA_Agraz3, lambda = seq(-10, 20, 1/10))
# lambda.MDA3 <- b$x[which.max(b$y)]
# Transformacion seleccionada
lambda.MDA3 <- 0.1
### Anova MDA
# Creación del modelo ANOVA de dos factores tipo III
mod_MDA3 <- aov( (M2.2)^lambda.MDA3 ~ N * P, data = MDA_Agraz3)
AOV_MDA3 <- Anova(mod_MDA3, type = "III"); AOV_MDA3## Anova Table (Type III tests)
##
## Response: (M2.2)^lambda.MDA3
## Sum Sq Df F value Pr(>F)
## (Intercept) 6.6598 1 1114.1845 < 2e-16 ***
## N 0.0696 2 5.8206 0.01123 *
## P 0.0002 2 0.0145 0.98560
## N:P 0.0079 4 0.3312 0.85336
## Residuals 0.1076 18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Residuales modelo bifactorial
res_MDA3 <- mod_MDA3$residuals
res.est_MDA3 <-res_MDA3/sd(res_MDA3)
# Probamos los supuestos
par(mfrow=c(2,2))
hist( (MDA_Agraz3$M2.2), main="Variable original" )
hist( (MDA_Agraz3$M2.2)^lambda.MDA3, main="Variable transformada" )
boxplot( (MDA_Agraz3$M2.2), horizontal = T, main="Variable original" )
boxplot( (MDA_Agraz3$M2.2)^lambda.MDA3, horizontal = T, main="Variable transformada" )
# Prueba de normalidad de los residuos
par(mfrow=c(1,2))
plot(mod_MDA3, which = 2)
library(nortest)
ad.test(res.est_MDA3)##
## Anderson-Darling normality test
##
## data: res.est_MDA3
## A = 0.31497, p-value = 0.5236cvm.test(res.est_MDA3)##
## Cramer-von Mises normality test
##
## data: res.est_MDA3
## W = 0.05021, p-value = 0.4953pearson.test(res.est_MDA3)##
## Pearson chi-square normality test
##
## data: res.est_MDA3
## P = 8.8519, p-value = 0.1151ks.test(res.est_MDA3, "pnorm")##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: res.est_MDA3
## D = 0.11951, p-value = 0.7925
## alternative hypothesis: two-sidedshapiro.test(res.est_MDA3)##
## Shapiro-Wilk normality test
##
## data: res.est_MDA3
## W = 0.96916, p-value = 0.5798# Prueba de homogeneidad de las varianzas (homocedasticidad)
plot(mod_MDA3, which = 1)
# Test con seguridad de normalidad (Bartlett, F-test) https://rpubs.com/Joaquin_AR/218466
bartlett.test( (M2.2)^lambda.MDA3 ~ N, data = MDA_Agraz3)##
## Bartlett test of homogeneity of variances
##
## data: (M2.2)^lambda.MDA3 by N
## Bartlett's K-squared = 2.9437, df = 2, p-value = 0.2295bartlett.test( (M2.2)^lambda.MDA3 ~ P, data = MDA_Agraz3)##
## Bartlett test of homogeneity of variances
##
## data: (M2.2)^lambda.MDA3 by P
## Bartlett's K-squared = 0.45062, df = 2, p-value = 0.7983# Analisis parametrico (Tukey 06 DDT)
library(agricolae)
Tuk_MD3 <- HSD.test(mod_MDA3, "N", console=F, alpha = 0.05)
Tukey_MDA3 <- cbind(
rownames_to_column(data.frame( Tukey=Tuk_MD3$groups[2] ), var = "N")
)# Grafico MDA 3
Resumen_MDA3 <- MDA_Agraz3 %>%
dplyr::group_by(., N) %>%
dplyr::summarise( n=n(), media=mean(M2.2), sd=sd(M2.2)) %>%
dplyr::mutate( se=sd/sqrt(n)) %>%
dplyr::mutate( ic=se * qt((1-0.05)/2 + .5, n-1)) %>%
dplyr::full_join(., Tukey_MDA3, by="N")
head(Resumen_MDA3)## # A tibble: 3 × 7
## N n media sd se ic groups
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 0.00 N 9 62.1 32.5 10.8 25.0 a
## 2 1.10 N 9 17.0 10.1 3.37 7.77 b
## 3 1.80 N 9 19.0 6.23 2.08 4.79 bFigura_MDA3 <- ggplot(Resumen_MDA3, aes(x=N, y=media) )+
geom_errorbar( aes(x=N, ymin=media-se, ymax=media+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=N, y=media)) +
geom_text( aes(label=groups, x=N, y=media+se+5 ), color="red", size=4)+
labs(x = "Niveles de nitrogeno",
y = (expression("MDA (mmol"~" g"^{-1}~")"))) +
ylim(0.0, 100)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_classic()
Figura_MDA3
### Kruskall MDA
# Analisis no parametrico
library(agricolae) ## Notas interesantes sobre Kruskall https://rpubs.com/Joaquin_AR/219504
MDA3.TTO <- kruskal.test(MDA_Agraz3$M2.2, MDA_Agraz3$TTO); MDA3.TTO##
## Kruskal-Wallis rank sum test
##
## data: MDA_Agraz3$M2.2 and MDA_Agraz3$TTO
## Kruskal-Wallis chi-squared = 17.333, df = 8, p-value = 0.02682MDA3.N <- kruskal.test(MDA_Agraz3$M2.2, MDA_Agraz3$N); MDA3.N##
## Kruskal-Wallis rank sum test
##
## data: MDA_Agraz3$M2.2 and MDA_Agraz3$N
## Kruskal-Wallis chi-squared = 16.152, df = 2, p-value = 0.000311MDA3.P <- kruskal.test(MDA_Agraz3$M2.2, MDA_Agraz3$TTO); MDA3.P##
## Kruskal-Wallis rank sum test
##
## data: MDA_Agraz3$M2.2 and MDA_Agraz3$TTO
## Kruskal-Wallis chi-squared = 17.333, df = 8, p-value = 0.02682# ver comparaciones puntuales entre tratamientos
KWlt.MDA3.TTO <- kruskal(MDA_Agraz3$M2.2, MDA_Agraz3$TTO, console=T,
p.adj="bonferroni",group=T)##
## Study: MDA_Agraz3$M2.2 ~ MDA_Agraz3$TTO
## Kruskal-Wallis test's
## Ties or no Ties
##
## Critical Value: 17.33333
## Degrees of freedom: 8
## Pvalue Chisq : 0.02681933
##
## MDA_Agraz3$TTO, means of the ranks
##
## MDA_Agraz3.M2.2 r
## 0.00 N/0.00 P 22.000000 3
## 0.00 N/0.25 P 22.333333 3
## 0.00 N/0.50 P 23.666667 3
## 1.10 N/0.00 P 8.000000 3
## 1.10 N/0.25 P 11.000000 3
## 1.10 N/0.50 P 8.666667 3
## 1.80 N/0.00 P 12.333333 3
## 1.80 N/0.25 P 11.333333 3
## 1.80 N/0.50 P 6.666667 3
##
## Post Hoc Analysis
##
## P value adjustment method: bonferroni
## t-Student: 3.774234
## Alpha : 0.05
## Minimum Significant Difference: 16.9724
##
## Treatments with the same letter are not significantly different.
##
## MDA_Agraz3$M2.2 groups
## 0.00 N/0.50 P 23.666667 a
## 0.00 N/0.25 P 22.333333 ab
## 0.00 N/0.00 P 22.000000 ab
## 1.80 N/0.00 P 12.333333 ab
## 1.80 N/0.25 P 11.333333 ab
## 1.10 N/0.25 P 11.000000 ab
## 1.10 N/0.50 P 8.666667 ab
## 1.10 N/0.00 P 8.000000 ab
## 1.80 N/0.50 P 6.666667 bKWlt.MDA3.N <- kruskal(MDA_Agraz3$M2.2, MDA_Agraz3$N, console=T,
p.adj="bonferroni",group=T)##
## Study: MDA_Agraz3$M2.2 ~ MDA_Agraz3$N
## Kruskal-Wallis test's
## Ties or no Ties
##
## Critical Value: 16.15168
## Degrees of freedom: 2
## Pvalue Chisq : 0.0003109627
##
## MDA_Agraz3$N, means of the ranks
##
## MDA_Agraz3.M2.2 r
## 0.00 N 22.666667 9
## 1.10 N 9.222222 9
## 1.80 N 10.111111 9
##
## Post Hoc Analysis
##
## P value adjustment method: bonferroni
## t-Student: 2.573641
## Alpha : 0.05
## Minimum Significant Difference: 6.168613
##
## Treatments with the same letter are not significantly different.
##
## MDA_Agraz3$M2.2 groups
## 0.00 N 22.666667 a
## 1.80 N 10.111111 b
## 1.10 N 9.222222 bKWlt.MDA3.P <- kruskal(MDA_Agraz3$M2.2, MDA_Agraz3$P, console=T,
p.adj="bonferroni",group=T)##
## Study: MDA_Agraz3$M2.2 ~ MDA_Agraz3$P
## Kruskal-Wallis test's
## Ties or no Ties
##
## Critical Value: 0.2574956
## Degrees of freedom: 2
## Pvalue Chisq : 0.8791957
##
## MDA_Agraz3$P, means of the ranks
##
## MDA_Agraz3.M2.2 r
## 0.00 P 14.11111 9
## 0.25 P 14.88889 9
## 0.50 P 13.00000 9
##
## Post Hoc Analysis
##
## P value adjustment method: bonferroni
## t-Student: 2.573641
## Alpha : 0.05
## Minimum Significant Difference: 9.973137
##
## Treatments with the same letter are not significantly different.
##
## MDA_Agraz3$M2.2 groups
## 0.25 P 14.88889 a
## 0.00 P 14.11111 a
## 0.50 P 13.00000 aGraficos jurados
#
Res_CCI2 <- rbind(Resumen_CCI2_07, Resumen_CCI2_16, Resumen_CCI2_25)
Res_CCI2$SDT <- as.factor(Res_CCI2$SDT)
ggplot(Res_CCI2, aes(x=TTO, y=median ) )+
facet_wrap(~SDT)+
geom_errorbar( aes(x=TTO, ymin=median-se, ymax=median+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=median) ) +
geom_text( aes(label=groups.kru, x=TTO, y=median+se+2.2 ), color="red", size=3.0)+
labs(x = "Tratamientos",
y = "Indice de contenido de clorofila (CCI)")+
ylim(0, 60)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_bw()+
theme(axis.text = element_text(angle = 90))+
theme(panel.grid = element_blank())
#
Res_CCI3 <- rbind(Resumen_CCI3_06, Resumen_CCI3_16, Resumen_CCI3_24)
Res_CCI3$SDT <- as.factor(Res_CCI3$SDT)
ggplot(Res_CCI3, aes(x=TTO, y=median ) )+
facet_wrap(~SDT)+
geom_errorbar( aes(x=TTO, ymin=median-se, ymax=median+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=median) ) +
geom_text( aes(label=groups.kru, x=TTO, y=median+se+2.2 ), color="red", size=3.0)+
labs(x = "Tratamientos",
y = "Indice de contenido de clorofila (CCI)")+
ylim(0, 60)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_bw()+
theme(axis.text = element_text(angle = 90))+
theme(panel.grid = element_blank())
#
Res_FC2 <- rbind(Resumen_FC2_07, Resumen_FC2_16, Resumen_FC2_25)
Res_FC2$SDT <- as.factor(Res_FC2$SDT)
ggplot(Res_FC2, aes(x=TTO, y=median ) )+
facet_wrap(~SDT)+
geom_errorbar( aes(x=TTO, ymin=median-se, ymax=median+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=median) ) +
geom_text( aes(label=groups.kru, x=TTO, y=median+se+0.03 ), color="red", size=3.0)+
labs(x = "Tratamientos",
y = "Eficiencia cuántica \nmáxima del PSII (Fv/Fm)")+
ylim(0.5, 1)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_bw()+
theme(axis.text = element_text(angle = 90))+
theme(panel.grid = element_blank())
#
Res_FC3 <- rbind(Resumen_FC3_06, Resumen_FC3_16, Resumen_FC3_24)
Res_FC3$SDT <- as.factor(Res_FC3$SDT)
ggplot(Res_FC3, aes(x=TTO, y=median ) )+
facet_wrap(~SDT)+
geom_errorbar( aes(x=TTO, ymin=median-se, ymax=median+se),
width=0.2, colour="gray", alpha=0.9)+
geom_point( aes(x=TTO, y=median) ) +
geom_text( aes(label=groups.kru, x=TTO, y=median+se+0.03 ), color="red", size=3.0)+
labs(x = "Tratamientos",
y = "Eficiencia cuántica \nmáxima del PSII (Fv/Fm)")+
ylim(0.5, 1)+
theme( axis.text.x=element_text(angle = 45, hjust = 2) )+
theme_bw()+
theme(axis.text = element_text(angle = 90))+
theme(panel.grid = element_blank())