setwd("~/Dropbox/Fernando 2018/archives/data")
library("lme4", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
## Loading required package: Matrix
enemy<- read.table("invas_fam_out_2.csv", sep = ",", header = T)
attach(enemy)
damage<-as.factor(damage)
fam<-as.factor(fam)
##Modelo sin interacciónes
fit1<-glmer(seeds.2~damage*origen+(1|pop), family = "poisson")
summary(fit1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: seeds.2 ~ damage * origen + (1 | pop)
##
## AIC BIC logLik deviance df.resid
## 5903.4 5923.0 -2946.7 5893.4 371
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.9370 -1.4710 0.1727 1.6322 15.7326
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop (Intercept) 2.629e-18 1.621e-09
## Number of obs: 376, groups: pop, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.78316 0.01857 203.764 < 2e-16 ***
## damage1 -0.15408 0.02819 -5.466 4.6e-08 ***
## origenMEX 0.27190 0.02206 12.324 < 2e-16 ***
## damage1:origenMEX 0.08254 0.03288 2.510 0.0121 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) damag1 orgMEX
## damage1 -0.659
## origenMEX -0.841 0.554
## dmg1:rgnMEX 0.565 -0.857 -0.671
anova(fit1)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## damage 1 29.2 29.2 29.1971
## origen 1 357.0 357.0 356.9980
## damage:origen 1 6.3 6.3 6.3019
##Probando efectos aleatorios
modsinpop<-glm(seeds.2~damage*origen, family="poisson")
anova(fit1, modsinpop)
## Data: NULL
## Models:
## modsinpop: seeds.2 ~ damage * origen
## fit1: seeds.2 ~ damage * origen + (1 | pop)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinpop 4 5901.4 5917.1 -2946.7 5893.4
## fit1 5 5903.4 5923.0 -2946.7 5893.4 0 1 1
#Finalmente el diagnostico de residuales del modelo...
plot(residuals(fit1)~predict(fit1,type="link"), xlab=expression(hat(eta)), ylab="Deviance residuals")
#...y los intervalos de confianza para las familias y poblaciones
library(lattice)
ranef(fit1)
## $pop
## (Intercept)
## morelos 4.763355e-17
## teotihuacan -4.763355e-17
## valdeflores -4.195194e-17
## zubia 4.195194e-17
dput(cl<-ranef(fit1,condVar=T))
## structure(list(pop = structure(list(`(Intercept)` = c(4.76335527703966e-17,
## -4.76335529799826e-17, -4.19519437746935e-17, 4.19519410692534e-17
## )), class = "data.frame", row.names = c("morelos", "teotihuacan",
## "valdeflores", "zubia"), postVar = structure(c(2.62921254120833e-18,
## 2.62921254120832e-18, 2.62921254120835e-18, 2.62921254120836e-18
## ), .Dim = c(1L, 1L, 4L)))), class = "ranef.mer")
dotplot(cl)
## $pop
