###Modelos Lineales Generalizados por Grupo Funcional
##Detrivoros
setwd("~/Dropbox/Proyecto_Neusa/data")
abun<-read.table("detrivo.csv",sep = ",",header = T)
attach(abun)
Zona<-as.factor(Zona)
Profundidad<-as.factor(Profundidad)
Month<-as.factor(Month)
###Modelo Lineal generalizado incluyendo todos los factores
library(lme4)
## Loading required package: Matrix
fit1<-glmer(Abundance~Month*Zona+(1|Profundidad),family = "poisson")
summary(fit1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Abundance ~ Month * Zona + (1 | Profundidad)
##
## AIC BIC logLik deviance df.resid
## 18305.3 18458.1 -9123.6 18247.3 1409
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.874 -1.499 -0.846 0.153 106.466
##
## Random effects:
## Groups Name Variance Std.Dev.
## Profundidad (Intercept) 0.2489 0.4989
## Number of obs: 1438, groups: Profundidad, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.39428 0.35776 12.283 < 2e-16 ***
## Month8 -4.25120 0.33919 -12.533 < 2e-16 ***
## Month14 -4.05789 0.19545 -20.762 < 2e-16 ***
## Month18 -4.34714 0.41336 -10.517 < 2e-16 ***
## Month22 -4.28785 0.25750 -16.652 < 2e-16 ***
## Month34 -4.70158 0.28440 -16.531 < 2e-16 ***
## Month39 -4.54074 0.27433 -16.552 < 2e-16 ***
## Zona2 -2.02765 0.06920 -29.303 < 2e-16 ***
## Zona3 -2.02583 0.07044 -28.760 < 2e-16 ***
## Zona4 -3.64489 0.08332 -43.745 < 2e-16 ***
## Month8:Zona2 3.52582 0.34506 10.218 < 2e-16 ***
## Month14:Zona2 2.88875 0.20509 14.085 < 2e-16 ***
## Month18:Zona2 4.05846 0.41772 9.716 < 2e-16 ***
## Month22:Zona2 3.06953 0.27748 11.062 < 2e-16 ***
## Month34:Zona2 4.00338 0.28879 13.863 < 2e-16 ***
## Month39:Zona2 3.65919 0.27914 13.109 < 2e-16 ***
## Month8:Zona3 2.48180 0.36197 6.856 7.07e-12 ***
## Month14:Zona3 2.58439 0.20871 12.382 < 2e-16 ***
## Month18:Zona3 2.30273 0.43732 5.266 1.40e-07 ***
## Month22:Zona3 2.53971 0.31121 8.161 3.33e-16 ***
## Month34:Zona3 2.97140 0.30653 9.694 < 2e-16 ***
## Month39:Zona3 3.00302 0.28656 10.479 < 2e-16 ***
## Month8:Zona4 4.48566 0.34913 12.848 < 2e-16 ***
## Month14:Zona4 3.96727 0.21165 18.744 < 2e-16 ***
## Month18:Zona4 4.45160 0.42013 10.596 < 2e-16 ***
## Month22:Zona4 4.42561 0.27030 16.373 < 2e-16 ***
## Month34:Zona4 5.03479 0.29550 17.038 < 2e-16 ***
## Month39:Zona4 4.20901 0.28954 14.537 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 28 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
##Modelo nulo
fit2<-glm(Abundance~1,family = "poisson")
##Significancia modelo completo
anova(fit1,fit2)
## Data: NULL
## Models:
## fit2: Abundance ~ 1
## fit1: Abundance ~ Month * Zona + (1 | Profundidad)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 1 23022 23027 -11509.9 23020
## fit1 29 18305 18458 -9123.6 18247 4772.6 28 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia de la interacción
drop1(fit1,test = "Chisq")
## Single term deletions
##
## Model:
## Abundance ~ Month * Zona + (1 | Profundidad)
## Df AIC LRT Pr(Chi)
## <none> 18305
## Month:Zona 18 19848 1579 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia del Mes
fit3<-glmer(Abundance~Zona+(1|Profundidad),family = "poisson")
anova(fit1, fit3)
## Data: NULL
## Models:
## fit3: Abundance ~ Zona + (1 | Profundidad)
## fit1: Abundance ~ Month * Zona + (1 | Profundidad)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit3 5 21264 21290 -10627.0 21254
## fit1 29 18305 18458 -9123.6 18247 3006.8 24 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia de Zona
fit4<-glmer(Abundance~Month+(1|Profundidad),family = "poisson")
anova(fit1,fit4)
## Data: NULL
## Models:
## fit4: Abundance ~ Month + (1 | Profundidad)
## fit1: Abundance ~ Month * Zona + (1 | Profundidad)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit4 8 21063 21105 -10523.4 21047
## fit1 29 18305 18458 -9123.6 18247 2799.5 21 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia de efecto aleatorio Profundidad
modsinprof<-glm(Abundance~Month*Zona+Grupo.Funcional,family = "poisson")
anova(fit1,modsinprof)
## Data: NULL
## Models:
## modsinprof: Abundance ~ Month * Zona + Grupo.Funcional
## fit1: Abundance ~ Month * Zona + (1 | Profundidad)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinprof 28 19245 19392 -9594.4 19189
## fit1 29 18305 18458 -9123.6 18247 941.45 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Diagnostico de residuales del modelo...
plot(residuals(fit1)~predict(fit1,type="link"), xlab=expression(hat(eta)), ylab="Deviance residuals")
# los intervalos de confianza para las profundidades
library(lattice)
dput(cl<-ranef(fit1,condVar=T))
## structure(list(Profundidad = structure(list(`(Intercept)` = c(0.498405020644126,
## -0.497791782253084)), class = "data.frame", row.names = c("1",
## "2"), postVar = structure(c(0.000156030217886417, 0.00104819959753633
## ), .Dim = c(1L, 1L, 2L)))), class = "ranef.mer")
dotplot(cl)
## $Profundidad
##Graficando las abundancias en función del tiempo con intervalos de confianza al 95% (inferencia directa)
fit1fit<-fitted(fit1)
abun$abundance.detrivo<-fit1fit
library(ggplot2)
p<-ggplot(abun,aes(Month,abundance.detrivo, group=Zona, color=Zona))
p + geom_smooth(method = "lm", se=F)+theme(legend.position="top")+stat_smooth(method = "lm",level = 0.95)
##Pruebas aposteriori interacción
Month<-as.factor(Month)
library(multcompView)
library(lsmeans)
## The 'lsmeans' package is being deprecated.
## Users are encouraged to switch to 'emmeans'.
## See help('transition') for more information, including how
## to convert 'lsmeans' objects and scripts to work with 'emmeans'.
marginal = lsmeans(fit1,
~ Zona:Month)
sum<-cld(marginal,
alpha = 0.05,
Letters = letters, ### Use lower-case letters for .group
adjust = "tukey") ### Tukey-adjusted comparisons
sum
## Zona Month lsmean SE df asymp.LCL asymp.UCL .group
## 1 34 -0.30729632 0.4492273 NA -1.70730541 1.092713 a
## 1 39 -0.14645664 0.4428089 NA -1.52646274 1.233549 abc
## 1 18 0.04714773 0.5400365 NA -1.63586646 1.730162 abcdefghi
## 1 22 0.10642912 0.4327406 NA -1.24219932 1.455058 abcde
## 1 8 0.14308673 0.4856741 NA -1.37050808 1.656682 abcdefgh
## 3 18 0.32404363 0.3793355 NA -0.85814880 1.506236 abc
## 1 14 0.33639047 0.3988015 NA -0.90646743 1.579248 abcde
## 4 39 0.41766508 0.3606268 NA -0.70622209 1.541552 ab
## 3 8 0.59905312 0.3732081 NA -0.56404332 1.762150 abcdefg
## 3 22 0.62030480 0.3924815 NA -0.60285690 1.843466 abcdefgh
## 3 34 0.63826994 0.3695014 NA -0.51327460 1.789814 abcdefg
## 4 14 0.65877118 0.3579904 NA -0.45689968 1.774442 abcd
## 4 0 0.74939438 0.3574855 NA -0.36470280 1.863492 bcdef
## 3 39 0.83073082 0.3610915 NA -0.29460458 1.956066 cdefg
## 4 18 0.85385860 0.3565703 NA -0.25738645 1.965104 cdefg
## 4 22 0.88715234 0.3582070 NA -0.22919351 2.003498 defg
## 3 14 0.89494800 0.3590871 NA -0.22414059 2.014037 defgh
## 4 8 0.98385407 0.3581138 NA -0.13220136 2.099909 efgh
## 4 34 1.08259925 0.3574085 NA -0.03125806 2.196457 gh
## 2 22 1.14830444 0.3665404 NA 0.00598767 2.290621 fghi
## 2 14 1.19749108 0.3570416 NA 0.08477737 2.310205 h
## 2 39 1.48508417 0.3554475 NA 0.37733829 2.592830 ij
## 2 8 1.64125175 0.3570916 NA 0.52838218 2.754121 j
## 2 34 1.66842988 0.3551086 NA 0.56174014 2.775120 j
## 2 18 2.07795864 0.3566195 NA 0.96656021 3.189357 k
## 2 0 2.36663040 0.3544564 NA 1.26197317 3.471288 l
## 3 0 2.36845432 0.3546990 NA 1.26304122 3.473867 l
## 1 0 4.39428405 0.3577556 NA 3.27934511 5.509223 m
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 28 estimates
## P value adjustment: tukey method for comparing a family of 28 estimates
## significance level used: alpha = 0.05