###Modelos Lineales Generalizados por Grupo Funcional
##Fitofagos
setwd("~/Dropbox/Proyecto_Neusa/data")
abun<-read.table("fitofa.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
## 5975.9 6126.5 -2959.0 5917.9 1302
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2053 -0.8704 -0.5369 0.1734 20.2954
##
## Random effects:
## Groups Name Variance Std.Dev.
## Profundidad (Intercept) 2.646e-18 1.627e-09
## Number of obs: 1331, groups: Profundidad, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.28768 0.50098 0.574 0.56581
## Month8 -0.28768 1.11780 -0.257 0.79690
## Month14 1.61656 0.52180 3.098 0.00195 **
## Month18 -0.28768 0.86618 -0.332 0.73979
## Month22 -0.28768 0.86641 -0.332 0.73986
## Month34 -0.12063 0.57268 -0.211 0.83317
## Month39 0.20875 0.54280 0.385 0.70054
## Zona2 0.49690 0.51060 0.973 0.33047
## Zona3 0.77844 0.50566 1.539 0.12370
## Zona4 0.21716 0.50991 0.426 0.67020
## Month8:Zona2 0.29156 1.12618 0.259 0.79572
## Month14:Zona2 -1.89371 0.54057 -3.503 0.00046 ***
## Month18:Zona2 0.47556 0.87660 0.543 0.58747
## Month22:Zona2 0.68580 0.87662 0.782 0.43403
## Month34:Zona2 -0.06975 0.58796 -0.119 0.90557
## Month39:Zona2 -0.23224 0.55828 -0.416 0.67742
## Month8:Zona3 -0.35193 1.12608 -0.313 0.75464
## Month14:Zona3 -1.61741 0.53164 -3.042 0.00235 **
## Month18:Zona3 0.22008 0.87226 0.252 0.80080
## Month22:Zona3 -0.24960 0.87934 -0.284 0.77653
## Month34:Zona3 -0.32898 0.58440 -0.563 0.57347
## Month39:Zona3 -0.42982 0.55323 -0.777 0.43720
## Month8:Zona4 0.31347 1.12557 0.279 0.78063
## Month14:Zona4 -1.18312 0.53446 -2.214 0.02685 *
## Month18:Zona4 0.87109 0.87336 0.997 0.31857
## Month22:Zona4 0.63692 0.87547 0.728 0.46691
## Month34:Zona4 0.67612 0.58469 1.156 0.24753
## Month39:Zona4 0.02105 0.55839 0.038 0.96992
## ---
## 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 6116.7 6121.9 -3057.4 6114.7
## fit1 29 5975.9 6126.5 -2959.0 5917.9 196.79 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> 5975.9
## Month:Zona 18 6083.1 143.19 < 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 6121.3 6147.3 -3055.7 6111.3
## fit1 29 5975.9 6126.5 -2959.0 5917.9 193.36 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 6079.1 6120.7 -3031.6 6063.1
## fit1 29 5975.9 6126.5 -2959.0 5917.9 145.22 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 5973.9 6119.3 -2959 5917.9
## fit1 29 5975.9 6126.5 -2959 5917.9 0 1 1
##Diagnostico de residuales del modelo...
plot(residuals(fit1)~predict(fit1,type="link"), xlab=expression(hat(eta)), ylab="Deviance residuals")
## Warning in sqrt(object$devResid()): NaNs produced
# los intervalos de confianza para las profundidades
library(lattice)
dput(cl<-ranef(fit1,condVar=T))
## structure(list(Profundidad = structure(list(`(Intercept)` = c(1.69125384381633e-17,
## -1.69125389602343e-17)), class = "data.frame", row.names = c("1",
## "2"), postVar = structure(c(2.6459772445922e-18, 2.64597724459221e-18
## ), .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.fitofa<-fit1fit
library(ggplot2)
p<-ggplot(abun,aes(Month,abundance.fitofa, 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 22 8.804069e-14 0.70698275 NA -2.20329924 2.2032992 abcdef
## 1 8 4.467537e-13 0.99945123 NA -3.11477209 3.1147721 abcdef
## 1 18 1.025958e-10 0.70687502 NA -2.20296352 2.2029635 abcdef
## 1 34 1.670541e-01 0.27736182 NA -0.69733913 1.0314473 abcde
## 1 0 2.876821e-01 0.50098166 NA -1.27361842 1.8489826 abcdef
## 3 8 4.265185e-01 0.11785058 NA 0.05923926 0.7937978 a
## 1 39 4.964369e-01 0.20854713 NA -0.15349657 1.1463703 abcde
## 4 0 5.048376e-01 0.09491696 NA 0.20903055 0.8006446 a
## 2 14 5.074300e-01 0.10101532 NA 0.19261758 0.8222425 ab
## 3 22 5.288441e-01 0.13363077 NA 0.11238621 0.9453021 abcd
## 4 8 5.306283e-01 0.09166985 NA 0.24494078 0.8163157 ab
## 2 34 5.942072e-01 0.08944292 NA 0.31545997 0.8729545 abc
## 3 34 6.165140e-01 0.09407191 NA 0.32334051 0.9096874 abc
## 4 39 7.346469e-01 0.09016779 NA 0.45364058 1.0156532 abcde
## 2 39 7.610978e-01 0.08543632 NA 0.49483707 1.0273586 abcde
## 2 0 7.845814e-01 0.09853493 NA 0.47749902 1.0916638 abcde
## 2 8 7.884574e-01 0.09534618 NA 0.49131268 1.0856020 abcde
## 3 39 8.450632e-01 0.08192307 NA 0.58975141 1.1003750 abcde
## 4 22 8.540775e-01 0.08219946 NA 0.59790438 1.1102507 abcde
## 4 14 9.382696e-01 0.06593806 NA 0.73277482 1.1437645 bcde
## 2 18 9.724610e-01 0.09167004 NA 0.68677295 1.2581491 abcde
## 3 18 9.985288e-01 0.07647180 NA 0.76020581 1.2368519 cde
## 4 34 1.060330e+00 0.06984306 NA 0.84266544 1.2779948 de
## 3 14 1.065276e+00 0.07516463 NA 0.83102663 1.2995251 de
## 3 0 1.066127e+00 0.06868009 NA 0.85208655 1.2801671 de
## 4 18 1.088250e+00 0.05892550 NA 0.90460923 1.2718898 e
## 2 22 1.182695e+00 0.08980250 NA 0.90282751 1.4625633 e
## 1 14 1.904237e+00 0.14586637 NA 1.44964750 2.3588274 f
##
## 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