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+(1|pop)+(1|pop:fam), family = "poisson")
summary(fit1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: seeds.2 ~ damage + (1 | pop) + (1 | pop:fam)
##
## AIC BIC logLik deviance df.resid
## 5316.5 5332.2 -2654.3 5308.5 372
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.0189 -1.2915 0.1637 1.3344 11.8492
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 0.04283 0.2070
## pop (Intercept) 0.01759 0.1326
## Number of obs: 376, groups: pop:fam, 52; pop, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.90547 0.07415 52.672 < 2e-16 ***
## damage1 -0.10053 0.01470 -6.839 7.98e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## damage1 -0.093
##Probando efectos fijos (damage) con la correción Kenward-Roger para datos desbalanceados
library(pbkrtest)
modtodos<-lmer(seeds.2~damage+(1|pop)+(1|pop:fam), REML = F)
modsindamage<-lmer(seeds.2~1+(1|pop)+(1|pop:fam), REML = F)
p.damage<-KRmodcomp(modtodos,modsindamage)
p.damage
## F-test with Kenward-Roger approximation; computing time: 0.15 sec.
## large : seeds.2 ~ damage + (1 | pop) + (1 | pop:fam)
## small : seeds.2 ~ 1 + (1 | pop) + (1 | pop:fam)
## stat ndf ddf F.scaling p.value
## Ftest 5.1432 1.0000 331.8785 1 0.02398 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Podemos mejorar la presición de la estimación con un boostrap paramétrico
pboost.damage<-PBmodcomp(modtodos,modsindamage) #Demora un minuto o dos dependiendo del equipo
summary (pboost.damage)
## Parametric bootstrap test; time: 24.97 sec; samples: 1000 extremes: 28;
## large : seeds.2 ~ damage + (1 | pop) + (1 | pop:fam)
## small : seeds.2 ~ 1 + (1 | pop) + (1 | pop:fam)
## stat df ddf p.value
## PBtest 5.0879 0.02897 *
## Gamma 5.0879 0.03280 *
## Bartlett 4.7359 1.0000 0.02954 *
## F 5.0879 1.0000 28.906 0.03184 *
## LRT 5.0879 1.0000 0.02409 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Probando efectos aleatorios
mod.sin.pop<-glmer(seeds.2~damage+(1|pop:fam), family="poisson")
mod.sin.fam<-glmer(seeds.2~damage+(1|pop), family="poisson")
#Se calculan las probabilidades asociadas a los efectos aleatorios, son las mas conservadoras sin usar correcciones (a la baja) del valor de p, debido a la naturaleza desbalanceada de los datos.
anova(fit1,mod.sin.pop)
## Data: NULL
## Models:
## mod.sin.pop: seeds.2 ~ damage + (1 | pop:fam)
## fit1: seeds.2 ~ damage + (1 | pop) + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.sin.pop 3 5322.6 5334.4 -2658.3 5316.6
## fit1 4 5316.5 5332.2 -2654.3 5308.5 8.1102 1 0.004402 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fit1,mod.sin.fam)
## Data: NULL
## Models:
## mod.sin.fam: seeds.2 ~ damage + (1 | pop)
## fit1: seeds.2 ~ damage + (1 | pop) + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.sin.fam 3 5927.8 5939.6 -2960.9 5921.8
## fit1 4 5316.5 5332.2 -2654.3 5308.5 613.25 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 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)
dotplot(ranef(fit1,condVar=TRUE))
## $`pop:fam`
##
## $pop
