###Modelos Lineales Generalizados por Grupo Funcional
##Predadores
setwd("~/Dropbox/Proyecto_Neusa/data")
abun<-read.table("preda.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")
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00135867
## (tol = 0.001, component 1)
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
## 6334.6 6496.0 -3138.3 6276.6 1905
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.8678 -0.6071 -0.4183 0.1220 11.5798
##
## Random effects:
## Groups Name Variance Std.Dev.
## Profundidad (Intercept) 0.01629 0.1276
## Number of obs: 1934, groups: Profundidad, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.514880 0.241132 2.135 0.0327 *
## Month8 -0.210652 0.311900 -0.675 0.4994
## Month14 0.034688 0.261771 0.133 0.8946
## Month18 -0.582456 0.417268 -1.396 0.1628
## Month22 -0.087236 0.305548 -0.286 0.7753
## Month34 -0.282200 0.294173 -0.959 0.3374
## Month39 -0.410470 0.319783 -1.284 0.1993
## Zona2 -0.017332 0.242665 -0.071 0.9431
## Zona3 -0.235751 0.247230 -0.954 0.3403
## Zona4 0.118702 0.230010 0.516 0.6058
## Month8:Zona2 0.364701 0.336320 1.084 0.2782
## Month14:Zona2 -0.209054 0.295872 -0.707 0.4798
## Month18:Zona2 0.476084 0.437892 1.087 0.2769
## Month22:Zona2 0.073125 0.338799 0.216 0.8291
## Month34:Zona2 0.503854 0.316436 1.592 0.1113
## Month39:Zona2 0.458719 0.340428 1.347 0.1778
## Month8:Zona3 0.288992 0.351342 0.823 0.4108
## Month14:Zona3 -0.098592 0.308398 -0.320 0.7492
## Month18:Zona3 0.595120 0.448290 1.328 0.1843
## Month22:Zona3 0.016807 0.357459 0.047 0.9625
## Month34:Zona3 0.189336 0.350646 0.540 0.5892
## Month39:Zona3 0.606032 0.352490 1.719 0.0856 .
## Month8:Zona4 0.003146 0.323024 0.010 0.9922
## Month14:Zona4 -0.220765 0.276289 -0.799 0.4243
## Month18:Zona4 0.578669 0.424199 1.364 0.1725
## Month22:Zona4 -0.029436 0.316861 -0.093 0.9260
## Month34:Zona4 -0.050732 0.310368 -0.163 0.8702
## Month39:Zona4 0.234463 0.335444 0.699 0.4846
## ---
## 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
## convergence code: 0
## Model failed to converge with max|grad| = 0.00135867 (tol = 0.001, component 1)
##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 6372.9 6378.5 -3185.5 6370.9
## fit1 29 6334.6 6496.0 -3138.3 6276.6 94.318 28 4.138e-09 ***
## ---
## 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> 6334.6
## Month:Zona 18 6340.5 41.924 0.001133 **
## ---
## 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 6334.3 6362.1 -3162.1 6324.3
## fit1 29 6334.6 6496.0 -3138.3 6276.6 47.665 24 0.002776 **
## ---
## 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 6357.6 6402.1 -3170.8 6341.6
## fit1 29 6334.6 6496.0 -3138.3 6276.6 65.014 21 2.163e-06 ***
## ---
## 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 6354.3 6510.1 -3149.1 6298.3
## fit1 29 6334.6 6496.0 -3138.3 6276.6 21.678 1 3.224e-06 ***
## ---
## 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.123477440850494,
## -0.122535586320735)), class = "data.frame", row.names = c("1",
## "2"), postVar = structure(c(0.000336493739160976, 0.00181511229377092
## ), .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.preda<-fit1fit
library(ggplot2)
p<-ggplot(abun,aes(Month,abundance.preda, 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 18 -0.06757673 0.3647885 NA -1.20443377 1.0692803 ab
## 1 39 0.10440998 0.2469328 NA -0.66515179 0.8739718 ab
## 3 34 0.18626437 0.1832988 NA -0.38498306 0.7575118 ab
## 3 22 0.20870051 0.1772068 NA -0.34356124 0.7609623 ab
## 3 14 0.21522442 0.1536693 NA -0.26368319 0.6941320 ab
## 1 34 0.23267970 0.2129944 NA -0.43111369 0.8964731 ab
## 3 0 0.27912862 0.1411055 NA -0.16062403 0.7188813 ab
## 3 18 0.29179244 0.1548097 NA -0.19066937 0.7742542 ab
## 4 34 0.30064941 0.1220837 NA -0.07982240 0.6811212 a
## 1 8 0.30422791 0.2362262 NA -0.43196700 1.0404228 ab
## 2 14 0.32318186 0.1355585 NA -0.09928376 0.7456475 ab
## 3 8 0.35746878 0.1521089 NA -0.11657590 0.8315135 ab
## 2 18 0.39117537 0.1307453 NA -0.01628994 0.7986407 ab
## 4 8 0.42607650 0.1105998 NA 0.08139426 0.7707587 ab
## 1 22 0.42764414 0.2271700 NA -0.28032729 1.1356156 ab
## 4 14 0.44750521 0.1142393 NA 0.09148041 0.8035300 ab
## 4 39 0.45757554 0.1241395 NA 0.07069686 0.8444542 ab
## 3 39 0.47469048 0.1384939 NA 0.04307659 0.9063044 ab
## 2 22 0.48343676 0.1440426 NA 0.03453066 0.9323429 ab
## 2 0 0.49754760 0.1333244 NA 0.08204439 0.9130508 ab
## 1 0 0.51487970 0.2411320 NA -0.23660404 1.2663634 ab
## 4 22 0.51691052 0.1103917 NA 0.17287689 0.8609442 ab
## 2 39 0.54579730 0.1139603 NA 0.19064198 0.9009526 ab
## 1 14 0.54956780 0.1647996 NA 0.03597272 1.0631629 ab
## 4 18 0.62979397 0.1053245 NA 0.30155207 0.9580359 ab
## 4 0 0.63358180 0.1077805 NA 0.29768565 0.9694779 ab
## 2 8 0.65159724 0.1226708 NA 0.26929580 1.0338987 ab
## 2 34 0.71920186 0.1137693 NA 0.36464188 1.0737618 b
##
## 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