setwd("/Volumes/GoogleDrive/Mi unidad/Agrosavia/Colaboraciones/Laura")
datos6<-read.table("sanluisflo.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
## -9.428 -5.720 -1.716 2.256 15.621
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.44113 0.04218 81.577 < 2e-16 ***
## forestalRoble -0.42614 0.06711 -6.350 2.16e-10 ***
## forestalTerminalia 0.09499 0.05829 1.630 0.103184
## genTCS 01 -0.93334 0.07940 -11.755 < 2e-16 ***
## genTCS 06 -0.38870 0.06636 -5.857 4.71e-09 ***
## genTCS13 -0.66508 0.07238 -9.188 < 2e-16 ***
## genTCS19 -0.06241 0.06061 -1.030 0.303168
## forestalRoble:genTCS 01 0.79317 0.11029 7.192 6.39e-13 ***
## forestalTerminalia:genTCS 01 0.64458 0.10041 6.420 1.37e-10 ***
## forestalRoble:genTCS 06 0.09281 0.10388 0.893 0.371629
## forestalTerminalia:genTCS 06 0.64683 0.08528 7.585 3.32e-14 ***
## forestalRoble:genTCS13 0.66508 0.10339 6.433 1.25e-10 ***
## forestalTerminalia:genTCS13 0.33666 0.09542 3.528 0.000418 ***
## forestalRoble:genTCS19 -0.03481 0.09696 -0.359 0.719553
## forestalTerminalia:genTCS19 -0.67711 0.09317 -7.267 3.67e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 8665.4 on 269 degrees of freedom
## Residual deviance: 7914.4 on 255 degrees of freedom
## AIC: 8861.2
##
## 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 269 8665.4
## forestal 2 226.23 267 8439.2 < 2.2e-16 ***
## gen 4 175.31 263 8263.9 < 2.2e-16 ***
## forestal:gen 8 349.44 255 7914.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
## -3.4157 -2.0276 -1.4907 0.4377 15.5293
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.6314 0.1043 15.648 < 2e-16 ***
## forestalRoble -0.6716 0.1793 -3.746 0.000180 ***
## forestalTerminalia 0.1322 0.1428 0.926 0.354686
## genTCS 01 -1.1896 0.2158 -5.512 3.56e-08 ***
## genTCS 06 -1.3029 0.2255 -5.777 7.61e-09 ***
## genTCS13 -1.1545 0.2130 -5.421 5.92e-08 ***
## genTCS19 -0.4443 0.1668 -2.664 0.007732 **
## forestalRoble:genTCS 01 1.6702 0.2846 5.868 4.42e-09 ***
## forestalTerminalia:genTCS 01 1.1408 0.2571 4.437 9.13e-06 ***
## forestalRoble:genTCS 06 0.4485 0.3495 1.283 0.199395
## forestalTerminalia:genTCS 06 0.2599 0.2957 0.879 0.379444
## forestalRoble:genTCS13 -0.2978 0.3969 -0.750 0.453138
## forestalTerminalia:genTCS13 -3.4995 1.0271 -3.407 0.000656 ***
## forestalRoble:genTCS19 0.4008 0.2671 1.501 0.133434
## forestalTerminalia:genTCS19 -0.5988 0.2537 -2.360 0.018269 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2072.0 on 269 degrees of freedom
## Residual deviance: 1744.5 on 255 degrees of freedom
## AIC: 2114.1
##
## 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 269 2072.0
## forestal 2 13.882 267 2058.1 0.0009671 ***
## gen 4 214.073 263 1844.0 < 2.2e-16 ***
## forestal:gen 8 99.524 255 1744.5 < 2.2e-16 ***
## ---
## 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
## CCN51 3.33 0.0261 Inf 3.28 3.38 A
## TCS 06 3.19 0.0289 Inf 3.13 3.25 B
## TCS19 3.03 0.0303 Inf 2.97 3.09 C
## TCS13 3.00 0.0306 Inf 2.94 3.06 C
## TCS 01 2.88 0.0330 Inf 2.81 2.94 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 3.32 0.0207 Inf 3.28 3.36 A
## Abarco 3.03 0.0239 Inf 2.98 3.08 B
## Roble 2.91 0.0247 Inf 2.86 2.96 C
##
## 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.79 0.0354 Inf 3.72 3.86 A
## Terminalia CCN51 3.54 0.0402 Inf 3.46 3.61 B
## Abarco CCN51 3.44 0.0422 Inf 3.36 3.52 BC
## Abarco TCS19 3.38 0.0435 Inf 3.29 3.46 BCD
## Terminalia TCS 01 3.25 0.0465 Inf 3.16 3.34 CDE
## Terminalia TCS13 3.21 0.0474 Inf 3.11 3.30 DE
## Abarco TCS 06 3.05 0.0512 Inf 2.95 3.15 EF
## Roble TCS13 3.01 0.0522 Inf 2.91 3.12 EFG
## Roble CCN51 3.01 0.0522 Inf 2.91 3.12 EFG
## Roble TCS19 2.92 0.0548 Inf 2.81 3.03 FGH
## Roble TCS 01 2.87 0.0560 Inf 2.77 2.98 FGH
## Terminalia TCS19 2.80 0.0582 Inf 2.68 2.91 FGHI
## Abarco TCS13 2.78 0.0588 Inf 2.66 2.89 GHI
## Roble TCS 06 2.72 0.0605 Inf 2.60 2.84 HI
## Abarco TCS 01 2.51 0.0673 Inf 2.38 2.64 I
##
## 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
## CCN51 1.452 0.0680 Inf 1.318 1.585 A
## TCS 01 1.199 0.0809 Inf 1.040 1.358 AB
## TCS19 0.941 0.0858 Inf 0.773 1.109 B
## TCS 06 0.385 0.1140 Inf 0.161 0.608 C
## TCS13 -0.969 0.3536 Inf -1.662 -0.276 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
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.813 0.0743 Inf 0.667542 0.959 A
## Roble 0.586 0.0889 Inf 0.411628 0.760 A
## Terminalia 0.406 0.2072 Inf -0.000336 0.812 A
##
## 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
## Terminalia CCN51 1.764 0.0976 Inf 1.5723 1.9549 A
## Terminalia TCS 01 1.715 0.1000 Inf 1.5188 1.9108 AB
## Abarco CCN51 1.631 0.1043 Inf 1.4271 1.8358 AB
## Roble TCS 01 1.440 0.1147 Inf 1.2155 1.6652 ABC
## Abarco TCS19 1.187 0.1302 Inf 0.9320 1.4423 BCD
## Roble CCN51 0.960 0.1459 Inf 0.6739 1.2457 CDE
## Roble TCS19 0.916 0.1491 Inf 0.6241 1.2085 CDE
## Terminalia TCS 06 0.721 0.1644 Inf 0.3983 1.0428 DE
## Terminalia TCS19 0.721 0.1644 Inf 0.3983 1.0428 DE
## Abarco TCS13 0.477 0.1857 Inf 0.1130 0.8409 DEF
## Abarco TCS 01 0.442 0.1890 Inf 0.0714 0.8122 DEF
## Abarco TCS 06 0.329 0.2000 Inf -0.0635 0.7205 EF
## Roble TCS 06 0.105 0.2236 Inf -0.3329 0.5436 EF
## Roble TCS13 -0.492 0.3015 Inf -1.0834 0.0985 F
## Terminalia TCS13 -2.890 1.0000 Inf -4.8503 -0.9304 F
##
## 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)