###Modelos Lineales Generalizados por Grupo Funcional
##Omnivoros
setwd("~/Dropbox/Proyecto_Neusa/data")
abun<-read.table("omni.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")
## fixed-effect model matrix is rank deficient so dropping 5 columns / coefficients
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
## 2416.6 2508.3 -1184.3 2368.6 314
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8201 -1.2383 -0.6829 0.0827 31.1838
##
## Random effects:
## Groups Name Variance Std.Dev.
## Profundidad (Intercept) 0.06239 0.2498
## Number of obs: 338, groups: Profundidad, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.42058 1.10293 0.381 0.70296
## Month8 -1.16709 0.61588 -1.895 0.05809 .
## Month14 -0.66412 1.47783 -0.449 0.65315
## Month18 -0.76131 0.39630 -1.921 0.05473 .
## Month22 0.43152 0.23568 1.831 0.06710 .
## Month34 -0.04549 1.23101 -0.037 0.97052
## Month39 -0.66393 0.43362 -1.531 0.12573
## Zona2 1.62363 1.09092 1.488 0.13667
## Zona3 0.07596 1.10093 0.069 0.94499
## Zona4 0.50284 1.06708 0.471 0.63747
## Month8:Zona2 0.91613 0.62744 1.460 0.14426
## Month14:Zona2 -0.39575 1.48450 -0.267 0.78979
## Month18:Zona2 -0.30818 0.42085 -0.732 0.46400
## Month22:Zona2 -2.23896 0.32028 -6.991 2.74e-12 ***
## Month34:Zona2 -0.93711 1.23673 -0.758 0.44861
## Month39:Zona2 -0.31783 0.45114 -0.705 0.48111
## Month8:Zona3 0.94519 0.66690 1.417 0.15640
## Month14:Zona3 1.27942 1.49196 0.858 0.39115
## Month18:Zona3 1.24213 0.44967 2.762 0.00574 **
## Month34:Zona3 -0.31889 1.27131 -0.251 0.80194
## Month39:Zona3 0.73959 0.48749 1.517 0.12923
## Month14:Zona4 -0.24061 1.51842 -0.158 0.87409
## Month34:Zona4 -0.83363 1.26389 -0.660 0.50953
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 23 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 5 columns / coefficients
##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 2713.1 2716.9 -1355.5 2711.1
## fit1 24 2416.6 2508.3 -1184.3 2368.6 342.51 23 < 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> 2416.6
## Month:Zona 13 2507.4 116.83 < 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 2586.1 2605.2 -1288.1 2576.1
## fit1 24 2416.6 2508.3 -1184.3 2368.6 207.56 19 < 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 2625.8 2656.4 -1304.9 2609.8
## fit1 24 2416.6 2508.3 -1184.3 2368.6 241.25 16 < 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 23 2443.7 2531.6 -1198.8 2397.7
## fit1 24 2416.6 2508.3 -1184.3 2368.6 29.128 1 6.773e-08 ***
## ---
## 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.243415613908621,
## -0.240320992761421)), class = "data.frame", row.names = c("1",
## "2"), postVar = structure(c(0.00106823110659227, 0.00571818527512103
## ), .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.omnivore<-fit1fit
library(ggplot2)
p<-ggplot(abun,aes(Month,abundance.omnivore, 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
## 4 8 -0.24366473 0.6054046 NA -2.13039741 1.643068 abcd
## 1 14 -0.24354374 1.0165966 NA -3.41174903 2.924662 abcde
## 1 39 -0.24335195 1.0146026 NA -3.40534295 2.918639 abcde
## 4 14 0.01868309 0.4458947 NA -1.37093974 1.408306 abc
## 4 34 0.04430761 0.3976861 NA -1.19507411 1.283689 abc
## 3 34 0.13216316 0.3226736 NA -0.87344334 1.137770 abc
## 4 18 0.16211384 0.3798440 NA -1.02166316 1.345891 abc
## 2 22 0.23676679 0.2685256 NA -0.60008842 1.073622 a
## 4 39 0.25948968 0.4173738 NA -1.04124850 1.560228 abc
## 3 8 0.27464329 0.2626267 NA -0.54382804 1.093115 ab
## 1 34 0.37509505 0.6032949 NA -1.50506295 2.255253 abcde
## 3 0 0.49654179 0.2457076 NA -0.26920170 1.262285 abc
## 3 39 0.57220345 0.2315267 NA -0.14934542 1.293752 abc
## 4 0 0.92342278 0.2775409 NA 0.05847145 1.788374 abc
## 3 22 0.92806562 0.2426061 NA 0.17198800 1.684143 abc
## 2 18 0.97471769 0.2093888 NA 0.32216129 1.627274 abc
## 3 18 0.97736022 0.2203096 NA 0.29076929 1.663951 abc
## 2 14 0.98433621 0.2131575 NA 0.32003464 1.648638 abc
## 2 34 1.06161647 0.1998912 NA 0.43865921 1.684574 bc
## 2 39 1.06243984 0.2023157 NA 0.43192639 1.692953 bc
## 3 14 1.11183420 0.2136120 NA 0.44611629 1.777552 c
## 2 8 1.79325534 0.1964839 NA 1.18091661 2.405594 de
## 2 0 2.04420793 0.1996024 NA 1.42215054 2.666265 e
## 1 0 nonEst NA NA NA NA
## 1 8 nonEst NA NA NA NA
## 1 18 nonEst NA NA NA NA
## 1 22 nonEst NA NA NA NA
## 4 22 nonEst NA NA NA NA
##
## 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