setwd("~/Dropbox/Proyecto_Neusa/data")
abun<-read.table("m10.csv",sep = ",",header = T)
attach(abun)
Zona<-as.factor(Zona)
Profundidad<-as.factor(Profundidad)
Month<-as.factor(Month)
Grupo.funcional<-as.factor(Grupo.Funcional)
library("lme4", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
## Loading required package: Matrix
###Modelo Lineal generalizado incluyendo todos los factores
fit1<-glmer(Abundance~Month*Zona+(1|Profundidad)+Grupo.Funcional,family = "poisson")
summary(fit1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Abundance ~ Month * Zona + (1 | Profundidad) + Grupo.Funcional
##
## AIC BIC logLik deviance df.resid
## 38197.2 38393.0 -19068.6 38137.2 5011
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.554 -1.051 -0.746 -0.061 185.546
##
## Random effects:
## Groups Name Variance Std.Dev.
## Profundidad (Intercept) 0.04876 0.2208
## Number of obs: 5041, groups: Profundidad, 2
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.130792 0.168764 18.551 < 2e-16 ***
## Month8 -2.434694 0.188517 -12.915 < 2e-16 ***
## Month14 -1.973392 0.104668 -18.854 < 2e-16 ***
## Month18 -2.793581 0.256356 -10.897 < 2e-16 ***
## Month22 -2.401543 0.166572 -14.417 < 2e-16 ***
## Month34 -2.708815 0.145356 -18.636 < 2e-16 ***
## Month39 -2.652768 0.144107 -18.408 < 2e-16 ***
## Zona2 -0.995574 0.063681 -15.634 < 2e-16 ***
## Zona3 -1.200836 0.064413 -18.643 < 2e-16 ***
## Zona4 -2.076339 0.067791 -30.628 < 2e-16 ***
## Grupo.Funcional -0.114499 0.008526 -13.429 < 2e-16 ***
## Month8:Zona2 1.848544 0.193967 9.530 < 2e-16 ***
## Month14:Zona2 1.029850 0.115303 8.932 < 2e-16 ***
## Month18:Zona2 2.195908 0.260569 8.427 < 2e-16 ***
## Month22:Zona2 1.428850 0.177489 8.050 8.26e-16 ***
## Month34:Zona2 2.099724 0.150825 13.922 < 2e-16 ***
## Month39:Zona2 1.934155 0.149813 12.910 < 2e-16 ***
## Month8:Zona3 1.243945 0.201601 6.170 6.81e-10 ***
## Month14:Zona3 1.267474 0.116614 10.869 < 2e-16 ***
## Month18:Zona3 1.862395 0.263496 7.068 1.57e-12 ***
## Month22:Zona3 1.284455 0.185552 6.922 4.44e-12 ***
## Month34:Zona3 1.608098 0.161027 9.986 < 2e-16 ***
## Month39:Zona3 1.758687 0.154033 11.418 < 2e-16 ***
## Month8:Zona4 2.369869 0.195880 12.099 < 2e-16 ***
## Month14:Zona4 1.935786 0.116688 16.589 < 2e-16 ***
## Month18:Zona4 2.917736 0.260673 11.193 < 2e-16 ***
## Month22:Zona4 2.441033 0.174666 13.975 < 2e-16 ***
## Month34:Zona4 2.840187 0.154651 18.365 < 2e-16 ***
## Month39:Zona4 2.486652 0.155740 15.967 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 29 > 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) + Grupo.Funcional
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 1 41381 41388 -20690 41379
## fit1 30 38197 38393 -19069 38137 3242.2 29 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia de la interacción y Grupo funcional
drop1(fit1,test = "Chisq")
## Single term deletions
##
## Model:
## Abundance ~ Month * Zona + (1 | Profundidad) + Grupo.Funcional
## Df AIC LRT Pr(Chi)
## <none> 38197
## Grupo.Funcional 1 38374 178.41 < 2.2e-16 ***
## Month:Zona 18 39165 1004.06 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia del Mes
fit3<-glmer(Abundance~Zona+(1|Profundidad)+Grupo.Funcional,family = "poisson")
anova(fit1, fit3)
## Data: NULL
## Models:
## fit3: Abundance ~ Zona + (1 | Profundidad) + Grupo.Funcional
## fit1: Abundance ~ Month * Zona + (1 | Profundidad) + Grupo.Funcional
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit3 6 40162 40201 -20075 40150
## fit1 30 38197 38393 -19069 38137 2012.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)+Grupo.Funcional,family = "poisson")
anova(fit1,fit4)
## Data: NULL
## Models:
## fit4: Abundance ~ Month + (1 | Profundidad) + Grupo.Funcional
## fit1: Abundance ~ Month * Zona + (1 | Profundidad) + Grupo.Funcional
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit4 9 39856 39914 -19919 39838
## fit1 30 38197 38393 -19069 38137 1700.4 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) + Grupo.Funcional
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinprof 29 38614 38803 -19278 38556
## fit1 30 38197 38393 -19069 38137 419.02 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.220278387217556,
## -0.220091370138298)), class = "data.frame", row.names = c("1",
## "2"), postVar = structure(c(8.01346763385051e-05, 0.000403547305662984
## ), .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<-fit1fit
library(ggplot2)
p<-ggplot(abun,aes(Month,abundance, group=Zona, color=Zona))
p + geom_smooth(method = "lm", se=F)+theme(legend.position="top")+stat_smooth(method = "lm",level = 0.95)
##Pruesbas 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.004626434 0.2947840 NA -0.92331556 0.9140627 abcdef
## 1 34 0.080139427 0.2055852 NA -0.56056316 0.7208420 a
## 1 39 0.136186565 0.2047086 NA -0.50178414 0.7741573 ab
## 1 8 0.354260728 0.2381779 NA -0.38801660 1.0965381 abcdef
## 1 22 0.387411667 0.2211263 NA -0.30172452 1.0765479 abcdef
## 3 8 0.397369940 0.1693054 NA -0.13026720 0.9250071 abc
## 3 22 0.471031115 0.1739718 NA -0.07114888 1.0132111 abcd
## 3 34 0.487401198 0.1685163 NA -0.03777675 1.0125791 abcd
## 4 39 0.546499864 0.1631108 NA 0.03816785 1.0548319 abcd
## 4 8 0.647790839 0.1610549 NA 0.14586600 1.1497157 bcde
## 3 18 0.656932765 0.1652254 NA 0.14201064 1.1718549 bcdef
## 4 14 0.675009443 0.1605720 NA 0.17458955 1.1754293 cdef
## 3 39 0.694037504 0.1630719 NA 0.18582696 1.2022480 cdef
## 4 0 0.712615702 0.1604035 NA 0.21272113 1.2125103 def
## 4 22 0.752106157 0.1608552 NA 0.25080368 1.2534086 def
## 1 14 0.815562501 0.1791472 NA 0.25725339 1.3738716 defg
## 2 22 0.820687634 0.1656140 NA 0.30455461 1.3368207 ef
## 4 18 0.836770547 0.1592354 NA 0.34051638 1.3330247 f
## 4 34 0.843987924 0.1608333 NA 0.34275383 1.3452220 ef
## 2 14 0.849838441 0.1613262 NA 0.34706807 1.3526088 ef
## 3 14 0.882200827 0.1619802 NA 0.37739224 1.3870094 fg
## 2 39 1.074767663 0.1592803 NA 0.57837328 1.5711620 gh
## 2 34 1.184289754 0.1590407 NA 0.68864225 1.6799373 h
## 2 18 1.195707534 0.1607568 NA 0.69471182 1.6967032 h
## 2 8 1.207231224 0.1604530 NA 0.70718230 1.7072801 h
## 3 0 1.588118529 0.1589322 NA 1.09280918 2.0834279 i
## 2 0 1.793380348 0.1586761 NA 1.29886921 2.2878915 j
## 1 0 2.788954282 0.1663396 NA 2.27055979 3.3073488 k
##
## 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