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 general
fit1<-glmer(seeds.2~damage*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 * pop + (1 | pop:fam)
##
## AIC BIC logLik deviance df.resid
## 5293.9 5329.2 -2637.9 5275.9 367
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.8497 -1.2407 0.1047 1.3524 11.2331
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 0.03942 0.1986
## Number of obs: 376, groups: pop:fam, 52
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.06637 0.05420 75.025 < 2e-16 ***
## damage1 -0.12609 0.02513 -5.018 5.23e-07 ***
## popteotihuacan -0.06870 0.07217 -0.952 0.341118
## popvaldeflores -0.31285 0.08656 -3.614 0.000301 ***
## popzubia -0.25374 0.09695 -2.617 0.008862 **
## damage1:popteotihuacan 0.09467 0.03446 2.747 0.006012 **
## damage1:popvaldeflores 0.01645 0.04450 0.370 0.711549
## damage1:popzubia -0.12818 0.05137 -2.495 0.012579 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) damag1 pptthc ppvldf popzub dmg1:ppt dmg1:ppv
## damage1 -0.225
## popteotihcn -0.751 0.169
## popvaldflrs -0.626 0.141 0.470
## popzubia -0.559 0.126 0.420 0.350
## dmg1:pptthc 0.164 -0.729 -0.239 -0.103 -0.092
## dmg1:ppvldf 0.127 -0.565 -0.095 -0.236 -0.071 0.412
## damag1:ppzb 0.110 -0.489 -0.083 -0.069 -0.226 0.357 0.276
##Modelo nulo
fit2<-glm(seeds.2~1, family="poisson")
##Significancia modelo completo
anova(fit1, fit2)
## Data: NULL
## Models:
## fit2: seeds.2 ~ 1
## fit1: seeds.2 ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 1 6309.2 6313.1 -3153.6 6307.2
## fit1 9 5293.9 5329.2 -2637.9 5275.9 1031.3 8 < 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:
## seeds.2 ~ damage * pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 5293.9
## damage:pop 3 5309.4 21.517 8.219e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia del daño
fit3<-glmer(seeds.2~pop+damage:pop+(1|pop:fam), family="poisson")
drop1(fit3, test = "Chisq")
## Single term deletions
##
## Model:
## seeds.2 ~ pop + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 5293.9
## pop:damage 4 5354.3 68.389 4.966e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia de la población
fit4<-glmer(seeds.2~damage+damage:pop+(1|pop:fam), family="poisson")
drop1(fit4, test = "Chisq")
## Single term deletions
##
## Model:
## seeds.2 ~ damage + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 5293.9
## damage:pop 6 5322.6 40.763 3.224e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia efecto aleatorio fam
modsinfam<-glm(seeds.2~damage*pop, family="poisson")
anova(fit1, modsinfam)
## Data: NULL
## Models:
## modsinfam: seeds.2 ~ damage * pop
## fit1: seeds.2 ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 8 5895.8 5927.2 -2939.9 5879.8
## fit1 9 5293.9 5329.2 -2637.9 5275.9 603.92 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Bootstrap paramétrico para confirmar significancia efecto fam (ho: la varianza es distinta de cero)
library(pbnm)
bootfam<-pbnm(fit1,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | pop:fam = 0
## from: glmer(formula = seeds.2 ~ damage * pop + (1 | pop:fam), family = "poisson")
## 1000 samples were taken Sun Oct 21 08:46:19 2018
## 2 samples had warnings, 2 in alternate model 0 in null model
## 2 unused samples. 0 <= P(abs((Intercept) | pop:fam) > |0.03942342|) <= 0.002
#Finalmente el diagnostico de residuales del modelo...
plot(residuals(fit1)~predict(fit1,type="link"), xlab=expression(hat(eta)), ylab="Deviance residuals")
##Graficando las normas de reacción para pop con intervalos de confianza al 95% (inferencia directa)
fit1fit<-fitted(fit1)
enemy$Seeds<-fit1fit
library(ggplot2)
p<-ggplot(enemy,aes(damage,Seeds, linetype=pop, color=pop))
p + geom_smooth(method = "lm", se=F)+scale_linetype_manual(values=c("twodash","solid","dotdash","dotted"))+theme(legend.position="top")+stat_smooth(method = "lm",level = 0.95)+scale_color_hue()
