setwd("/Volumes/GoogleDrive/Mi unidad/Agrosavia/Colaboraciones/Laura")
datos6<-read.table("santamariaflo.csv", header=T, sep=',')
datos6$gen<-as.factor(datos6$gen)
datos6$bloque<-as.factor(datos6$bloque)
datos6$forestal<-as.factor(datos6$forestal)
attach(datos6)
library(ggplot2)
library(Rmisc)
## Loading required package: lattice
## Loading required package: plyr
library(emmeans)
#Gráfica cojines por planta
ggplot(datos6, aes(fill=gen, y=cojplanta, x=forestal)) +
geom_bar(position="dodge", stat="identity") +
labs(x="Forestal", y="Cojines por planta")
# Gráfica frutos por planta
ggplot(datos6, aes(fill=gen, y=frutos, x=forestal)) +
geom_bar(position="dodge", stat="identity") +
labs(x="Forestal", y="Frutos por planta")
#modelos
gfit1<-glm(cojplanta~forestal*gen, data=datos6, family = poisson (link=log))
gfit2<-glm(frutos~forestal*gen, data=datos6, family = poisson (link=log))
# resultado cojines por planta
summary(gfit1)
##
## Call:
## glm(formula = cojplanta ~ forestal * gen, family = poisson(link = log),
## data = datos6)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.2269 -2.5166 -1.5811 0.8449 10.9016
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.49870 0.08276 30.192 < 2e-16 ***
## forestalRoble -1.34602 0.18211 -7.391 1.46e-13 ***
## forestalTerminalia 0.71349 0.10102 7.063 1.63e-12 ***
## genTCS 01 -2.09323 0.24981 -8.379 < 2e-16 ***
## genTCS 06 0.52578 0.10439 5.037 4.74e-07 ***
## genTCS13 -0.32965 0.12796 -2.576 0.00999 **
## genTCS19 -0.16333 0.12212 -1.337 0.18110
## forestalRoble:genTCS 01 1.45138 0.37245 3.897 9.74e-05 ***
## forestalTerminalia:genTCS 01 -0.55934 0.33671 -1.661 0.09668 .
## forestalRoble:genTCS 06 1.47854 0.20189 7.323 2.42e-13 ***
## forestalTerminalia:genTCS 06 0.18730 0.12609 1.485 0.13742
## forestalRoble:genTCS13 -0.59989 0.33069 -1.814 0.06967 .
## forestalTerminalia:genTCS13 -0.97300 0.17909 -5.433 5.54e-08 ***
## forestalRoble:genTCS19 1.47450 0.21982 6.708 1.98e-11 ***
## forestalTerminalia:genTCS19 -1.49072 0.18924 -7.877 3.35e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4213.2 on 179 degrees of freedom
## Residual deviance: 2173.4 on 165 degrees of freedom
## AIC: 2682.3
##
## Number of Fisher Scoring iterations: 6
anova(gfit1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: cojplanta
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 179 4213.2
## forestal 2 230.62 177 3982.5 < 2.2e-16 ***
## gen 4 1491.65 173 2490.9 < 2.2e-16 ***
## forestal:gen 8 317.54 165 2173.4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# resultado frutos por planta
summary(gfit2)
##
## Call:
## glm(formula = frutos ~ forestal * gen, family = poisson(link = log),
## data = datos6)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7795 -1.4287 -0.9129 0.0000 4.8796
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.877e-01 2.500e-01 1.151 0.2498
## forestalRoble -1.674e+00 6.292e-01 -2.661 0.0078 **
## forestalTerminalia -5.754e-01 4.167e-01 -1.381 0.1673
## genTCS 01 9.174e-17 3.536e-01 0.000 1.0000
## genTCS 06 -4.700e-01 4.031e-01 -1.166 0.2436
## genTCS13 -2.877e-01 3.819e-01 -0.753 0.4513
## genTCS19 1.719e-01 3.393e-01 0.506 0.6125
## forestalRoble:genTCS 01 1.386e+00 7.360e-01 1.884 0.0596 .
## forestalTerminalia:genTCS 01 -1.504e+00 8.580e-01 -1.753 0.0796 .
## forestalRoble:genTCS 06 1.936e+00 7.568e-01 2.559 0.0105 *
## forestalTerminalia:genTCS 06 9.118e-01 5.874e-01 1.552 0.1206
## forestalRoble:genTCS13 7.985e-01 8.241e-01 0.969 0.3326
## forestalTerminalia:genTCS13 3.637e-02 6.323e-01 0.058 0.9541
## forestalRoble:genTCS19 1.502e+00 7.148e-01 2.101 0.0356 *
## forestalTerminalia:genTCS19 -2.369e+00 1.107e+00 -2.139 0.0324 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 451.93 on 179 degrees of freedom
## Residual deviance: 399.68 on 165 degrees of freedom
## AIC: 574.88
##
## Number of Fisher Scoring iterations: 7
anova(gfit2, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: frutos
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 179 451.93
## forestal 2 15.683 177 436.25 0.0003931 ***
## gen 4 3.907 173 432.34 0.4187196
## forestal:gen 8 32.662 165 399.68 7.086e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Contrastes para Cojines
#Gen
contrast <- emmeans(gfit1, ~ gen)
## NOTE: Results may be misleading due to involvement in interactions
plot(contrast, comparisons = TRUE, xlab ="LN(Cojines por planta)")
cld_gen <-multcomp::cld(contrast, alpha = 0.05, Letters = LETTERS, reversed=T)
cld_gen
## gen emmean SE df asymp.LCL asymp.UCL .group
## TCS 06 3.369 0.0320 Inf 3.306 3.432 A
## CCN51 2.288 0.0637 Inf 2.163 2.413 B
## TCS19 2.119 0.0603 Inf 2.001 2.237 B
## TCS13 1.434 0.0992 Inf 1.240 1.628 C
## TCS 01 0.492 0.1305 Inf 0.236 0.748 D
##
## Results are averaged over the levels of: forestal
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 5 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
#Forestal
contrast <- emmeans(gfit1, ~ forestal)
## NOTE: Results may be misleading due to involvement in interactions
plot(contrast, comparisons = TRUE, xlab ="LN(Cojines por planta)")
cld_gen <-multcomp::cld(contrast, alpha = 0.05, Letters = LETTERS, reversed=T)
cld_gen
## forestal emmean SE df asymp.LCL asymp.UCL .group
## Terminalia 2.23 0.0574 Inf 2.12 2.35 A
## Abarco 2.09 0.0580 Inf 1.97 2.20 A
## Roble 1.50 0.0784 Inf 1.35 1.66 B
##
## Results are averaged over the levels of: gen
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
#Forestal*Gen
contrast <- emmeans(gfit1, ~ forestal * gen)
plot(contrast, comparisons = TRUE, xlab ="LN(Cojines por planta)")
cld_gen <-multcomp::cld(contrast, alpha = 0.05, Letters = LETTERS, reversed=T)
cld_gen
## forestal gen emmean SE df asymp.LCL asymp.UCL .group
## Terminalia TCS 06 3.925 0.0406 Inf 3.8458 4.005 A
## Terminalia CCN51 3.212 0.0579 Inf 3.0986 3.326 B
## Roble TCS 06 3.157 0.0595 Inf 3.0403 3.274 B
## Abarco TCS 06 3.024 0.0636 Inf 2.8998 3.149 B
## Abarco CCN51 2.499 0.0828 Inf 2.3365 2.661 C
## Roble TCS19 2.464 0.0842 Inf 2.2988 2.629 C
## Abarco TCS19 2.335 0.0898 Inf 2.1594 2.511 CD
## Abarco TCS13 2.169 0.0976 Inf 1.9778 2.360 CD
## Terminalia TCS13 1.910 0.1111 Inf 1.6918 2.127 DE
## Terminalia TCS19 1.558 0.1325 Inf 1.2985 1.818 EF
## Roble CCN51 1.153 0.1622 Inf 0.8347 1.471 FG
## Terminalia TCS 01 0.560 0.2182 Inf 0.1319 0.987 G
## Roble TCS 01 0.511 0.2236 Inf 0.0726 0.949 G
## Abarco TCS 01 0.405 0.2357 Inf -0.0565 0.867 G
## Roble TCS13 0.223 0.2582 Inf -0.2829 0.729 G
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 15 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
#Contrastes para frutos
#Gen
contrast2 <- emmeans(gfit2, ~ gen)
## NOTE: Results may be misleading due to involvement in interactions
plot(contrast2, comparisons = TRUE, xlab = "LN(Frutos por planta)")
cld_gen <-multcomp::cld(contrast2, alpha = 0.05, Letters = LETTERS, reversed=T)
cld_gen
## gen emmean SE df asymp.LCL asymp.UCL .group
## TCS 06 0.0173 0.166 Inf -0.308 0.34287 A
## CCN51 -0.4621 0.237 Inf -0.927 0.00307 A
## TCS13 -0.4715 0.218 Inf -0.898 -0.04498 A
## TCS 01 -0.5014 0.268 Inf -1.026 0.02367 A
## TCS19 -0.5792 0.352 Inf -1.269 0.11068 A
##
## Results are averaged over the levels of: forestal
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 5 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
#Forestal
contrast2 <- emmeans(gfit2, ~ forestal)
## NOTE: Results may be misleading due to involvement in interactions
plot(contrast2, comparisons = TRUE, xlab = "LN(Frutos por planta)")
cld_gen <-multcomp::cld(contrast2, alpha = 0.05, Letters = LETTERS, reversed=T)
cld_gen
## forestal emmean SE df asymp.LCL asymp.UCL .group
## Abarco 0.171 0.120 Inf -0.065 0.406 A
## Roble -0.379 0.174 Inf -0.720 -0.038 B
## Terminalia -0.990 0.270 Inf -1.519 -0.460 B
##
## Results are averaged over the levels of: gen
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
#Forestal*Gen
contrast2 <- emmeans(gfit2, ~ forestal * gen)
plot(contrast2, comparisons = TRUE, xlab = "LN(Frutos por planta)")
cld_gen <-multcomp::cld(contrast2, alpha = 0.05, Letters = LETTERS, reversed=T)
cld_gen
## forestal gen emmean SE df asymp.LCL asymp.UCL .group
## Abarco TCS19 0.460 0.229 Inf 0.00989 0.90918 A
## Abarco TCS 01 0.288 0.250 Inf -0.20231 0.77767 A
## Abarco CCN51 0.288 0.250 Inf -0.20231 0.77767 A
## Roble TCS19 0.288 0.250 Inf -0.20231 0.77767 A
## Terminalia TCS 06 0.154 0.267 Inf -0.36967 0.67797 A
## Roble TCS 06 0.080 0.277 Inf -0.46355 0.62364 A
## Roble TCS 01 0.000 0.289 Inf -0.56579 0.56579 A
## Abarco TCS13 0.000 0.289 Inf -0.56579 0.56579 A
## Abarco TCS 06 -0.182 0.316 Inf -0.80212 0.43747 A
## Terminalia CCN51 -0.288 0.333 Inf -0.94100 0.36564 A
## Terminalia TCS13 -0.539 0.378 Inf -1.27979 0.20180 A
## Roble TCS13 -0.875 0.447 Inf -1.75199 0.00105 A
## Roble CCN51 -1.386 0.577 Inf -2.51788 -0.25471 A
## Terminalia TCS 01 -1.792 0.707 Inf -3.17766 -0.40586 A
## Terminalia TCS19 -2.485 1.000 Inf -4.44487 -0.52494 A
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 15 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
detach(datos6)