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("todos_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
## 5026.9 5061.6 -2504.4 5008.9 342
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.9459 -1.2132 0.1518 1.4238 11.1978
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 0.0402 0.2005
## Number of obs: 351, groups: pop:fam, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.03226 0.05701 70.730 < 2e-16 ***
## damage1 -0.10931 0.02705 -4.041 5.33e-05 ***
## popteotihuacan -0.02886 0.07527 -0.383 0.70138
## popvaldeflores -0.27353 0.08885 -3.079 0.00208 **
## popzubia -0.21966 0.09911 -2.216 0.02667 *
## damage1:popteotihuacan 0.08099 0.03634 2.229 0.02583 *
## damage1:popvaldeflores 0.02171 0.04602 0.472 0.63711
## damage1:popzubia -0.14501 0.05233 -2.771 0.00559 **
## ---
## 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.242
## popteotihcn -0.757 0.183
## popvaldflrs -0.642 0.155 0.486
## popzubia -0.575 0.139 0.436 0.369
## dmg1:pptthc 0.180 -0.744 -0.243 -0.116 -0.104
## dmg1:ppvldf 0.142 -0.588 -0.108 -0.241 -0.082 0.438
## damag1:ppzb 0.125 -0.517 -0.095 -0.080 -0.230 0.385 0.304
##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 5905.3 5909.1 -2951.6 5903.3
## fit1 9 5026.9 5061.6 -2504.4 5008.9 894.4 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> 5026.9
## damage:pop 3 5041.4 20.5 0.0001337 ***
## ---
## 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> 5026.9
## pop:damage 4 5074.6 55.725 2.29e-11 ***
## ---
## 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> 5026.9
## damage:pop 6 5052.1 37.229 1.589e-06 ***
## ---
## 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 5581.0 5611.9 -2782.5 5565.0
## fit1 9 5026.9 5061.6 -2504.4 5008.9 556.16 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 Mon Oct 29 16:32:41 2018
## 1 samples had warnings, 1 in alternate model 0 in null model
## 1 unused samples. 0 <= P(abs((Intercept) | pop:fam) > |0.04020354|) <= 0.001
#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$Seeds2<-fit1fit
library(ggplot2)
p<-ggplot(enemy,aes(damage,Seeds2, 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)
##ESFUERZO REPRODUCTIVO: Modelo general
fit1<-lmer(arep~damage*pop+(1|pop:fam))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: arep ~ damage * pop + (1 | pop:fam)
##
## REML criterion at convergence: 2093.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1959 -0.5442 0.0933 0.6017 2.5686
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 6.764 2.601
## Residual 20.451 4.522
## Number of obs: 351, groups: pop:fam, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 33.4087 0.9573 34.900
## damage1 -0.3992 0.8993 -0.444
## popteotihuacan -1.2379 1.2669 -0.977
## popvaldeflores 1.8114 1.4622 1.239
## popzubia 0.7939 1.6273 0.488
## damage1:popteotihuacan 0.3644 1.2117 0.301
## damage1:popvaldeflores 0.3941 1.4053 0.280
## damage1:popzubia 0.6932 1.5673 0.442
##
## Correlation of Fixed Effects:
## (Intr) damag1 pptthc ppvldf popzub dmg1:ppt dmg1:ppv
## damage1 -0.503
## popteotihcn -0.756 0.380
## popvaldflrs -0.655 0.330 0.495
## popzubia -0.588 0.296 0.444 0.385
## dmg1:pptthc 0.374 -0.742 -0.498 -0.245 -0.220
## dmg1:ppvldf 0.322 -0.640 -0.243 -0.473 -0.189 0.475
## damag1:ppzb 0.289 -0.574 -0.218 -0.189 -0.473 0.426 0.367
##Modelo nulo
fit2<-lm(arep~1)
##Significancia modelo completo
anova(fit1, fit2)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## fit2: arep ~ 1
## fit1: arep ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 2163.4 2171.1 -1079.7 2159.4
## fit1 10 2126.9 2165.5 -1053.5 2106.9 52.464 8 1.369e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia interacción
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## arep ~ damage * pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 2126.9
## damage:pop 3 2121.1 0.22028 0.9743
##Significancia del daño
fit3<-lmer(arep~pop+damage:pop+(1|pop:fam))
drop1(fit3, test = "Chisq")
## Single term deletions
##
## Model:
## arep ~ pop + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 2126.9
## pop:damage 4 2119.2 0.25285 0.9927
##Significancia de la población
fit4<-lmer(arep~damage+damage:pop+(1|pop:fam))
drop1(fit4, test = "Chisq")
## Single term deletions
##
## Model:
## arep ~ damage + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 2126.9
## damage:pop 6 2122.3 7.3517 0.2895
##Significancia efecto aleatorio fam
modsinfam<-lm(arep~damage*pop)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: arep ~ damage * pop
## fit1: arep ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 9 2158.0 2192.8 -1070.0 2140.0
## fit1 10 2126.9 2165.5 -1053.5 2106.9 33.129 1 8.625e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Bootstrap paramétrico para reconfirmar significancia efecto del daño*fam (ho: la varianza es distinta de cero)
library(pbnm)
fit1f<-lmer(arep~damage*pop+(1|pop:fam), REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | pop:fam = 0
## from: lmer(formula = arep ~ damage * pop + (1 | pop:fam), REML = F)
## 1000 samples were taken Mon Oct 29 16:33:00 2018
## 2 samples had warnings, 2 in alternate model 0 in null model
## 2 unused samples. 0 <= P(abs((Intercept) | pop:fam) > |5.998127|) <= 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$arepfit2<-fit1fit
library(ggplot2)
p<-ggplot(enemy,aes(damage,arepfit2, 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)
##LMF: Modelo general
fit1<-lmer(almf~damage*pop+(1|pop:fam))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: almf ~ damage * pop + (1 | pop:fam)
##
## REML criterion at convergence: 1968
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4327 -0.5178 -0.0134 0.3997 8.0859
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 0.7146 0.8454
## Residual 16.0657 4.0082
## Number of obs: 351, groups: pop:fam, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 18.8663 0.6225 30.306
## damage1 -1.1648 0.7935 -1.468
## popteotihuacan -0.8608 0.8259 -1.042
## popvaldeflores -1.1429 0.9398 -1.216
## popzubia -1.5086 1.0425 -1.447
## damage1:popteotihuacan 0.8231 1.0689 0.770
## damage1:popvaldeflores 0.8914 1.2416 0.718
## damage1:popzubia 2.8967 1.3799 2.099
##
## Correlation of Fixed Effects:
## (Intr) damag1 pptthc ppvldf popzub dmg1:ppt dmg1:ppv
## damage1 -0.681
## popteotihcn -0.754 0.514
## popvaldflrs -0.662 0.451 0.499
## popzubia -0.597 0.407 0.450 0.396
## dmg1:pptthc 0.506 -0.742 -0.672 -0.335 -0.302
## dmg1:ppvldf 0.435 -0.639 -0.328 -0.651 -0.260 0.474
## damag1:ppzb 0.392 -0.575 -0.295 -0.260 -0.651 0.427 0.368
##Modelo nulo
fit2<-lm(almf~1)
##Significancia modelo completo
anova(fit1, fit2)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## fit2: almf ~ 1
## fit1: almf ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 1987.0 1994.8 -991.52 1983.0
## fit1 10 1996.1 2034.7 -988.03 1976.1 6.9795 8 0.5388
##Significancia interacción
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## almf ~ damage * pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 1996.1
## damage:pop 3 1994.5 4.4663 0.2153
##Significancia del daño
fit3<-lmer(almf~pop+damage:pop+(1|pop:fam))
drop1(fit3, test = "Chisq")
## Single term deletions
##
## Model:
## almf ~ pop + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 1996.1
## pop:damage 4 1992.9 4.8699 0.3009
##Significancia de la población
fit4<-lmer(almf~damage+damage:pop+(1|pop:fam))
drop1(fit4, test = "Chisq")
## Single term deletions
##
## Model:
## almf ~ damage + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 1996.1
## damage:pop 6 1989.7 5.6113 0.4681
##Significancia efecto aleatorio fam
modsinfam<-lm(almf~damage*pop)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: almf ~ damage * pop
## fit1: almf ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 9 1995.0 2029.7 -988.48 1977.0
## fit1 10 1996.1 2034.7 -988.03 1976.1 0.8977 1 0.3434
##Bootstrap paramétrico para reconfirmar significancia efecto del daño*fam (ho: la varianza es distinta de cero)
fit1f<-lmer(almf~damage*pop+(1|pop:fam), REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | pop:fam = 0
## from: lmer(formula = almf ~ damage * pop + (1 | pop:fam), REML = F)
## 1000 samples were taken Mon Oct 29 16:33:18 2018
## P(abs((Intercept) | pop:fam) > |0.4957279|) = 0.095
#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$almffit2<-fit1fit
p<-ggplot(enemy,aes(damage,almffit2, 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)
##SMF: Modelo general
fit1<-lmer(asmf~damage*pop+(1|pop:fam))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: asmf ~ damage * pop + (1 | pop:fam)
##
## REML criterion at convergence: 1970.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8634 -0.5732 -0.0079 0.4785 3.8842
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 4.416 2.101
## Residual 14.388 3.793
## Number of obs: 351, groups: pop:fam, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 34.0792 0.7874 43.280
## damage1 -0.3541 0.7542 -0.470
## popteotihuacan -0.2702 1.0423 -0.259
## popvaldeflores -0.6764 1.2021 -0.563
## popzubia -2.2380 1.3377 -1.673
## damage1:popteotihuacan 0.4336 1.0162 0.427
## damage1:popvaldeflores 0.8277 1.1786 0.702
## damage1:popzubia 0.5342 1.3143 0.406
##
## Correlation of Fixed Effects:
## (Intr) damag1 pptthc ppvldf popzub dmg1:ppt dmg1:ppv
## damage1 -0.513
## popteotihcn -0.755 0.388
## popvaldflrs -0.655 0.336 0.495
## popzubia -0.589 0.302 0.445 0.386
## dmg1:pptthc 0.381 -0.742 -0.507 -0.249 -0.224
## dmg1:ppvldf 0.328 -0.640 -0.248 -0.483 -0.193 0.475
## damag1:ppzb 0.294 -0.574 -0.222 -0.193 -0.483 0.426 0.367
##Modelo nulo
fit2<-lm(asmf~1)
##Significancia modelo completo
anova(fit1, fit2)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## fit2: asmf ~ 1
## fit1: asmf ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 2024.2 2032.0 -1010.12 2020.2
## fit1 10 2000.9 2039.5 -990.45 1980.9 39.332 8 4.264e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia interacción
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## asmf ~ damage * pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 2000.9
## damage:pop 3 1995.4 0.52217 0.914
##Significancia del daño
fit3<-lmer(asmf~pop+damage:pop+(1|pop:fam))
drop1(fit3, test = "Chisq")
## Single term deletions
##
## Model:
## asmf ~ pop + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 2000.9
## pop:damage 4 1993.4 0.53528 0.97
##Significancia de la población
fit4<-lmer(asmf~damage+damage:pop+(1|pop:fam))
drop1(fit4, test = "Chisq")
## Single term deletions
##
## Model:
## asmf ~ damage + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 2000.9
## damage:pop 6 1993.0 4.095 0.6638
##Significancia efecto aleatorio fam
modsinfam<-lm(asmf~damage*pop)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: asmf ~ damage * pop
## fit1: asmf ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 9 2029.2 2063.9 -1005.58 2011.2
## fit1 10 2000.9 2039.5 -990.45 1980.9 30.254 1 3.791e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Bootstrap paramétrico para reconfirmar significancia efecto del daño*fam (ho: la varianza es distinta de cero)
fit1f<-lmer(asmf~damage*pop+(1|pop:fam), REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | pop:fam = 0
## from: lmer(formula = asmf ~ damage * pop + (1 | pop:fam), REML = F)
## 1000 samples were taken Mon Oct 29 16:33:36 2018
## P(abs((Intercept) | pop:fam) > |3.906617|) = 0
#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$asmffit2<-fit1fit
p<-ggplot(enemy,aes(damage,asmffit2, 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)
##RMF: Modelo general
fit1<-lmer(armf~damage*pop+(1|pop:fam))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: armf ~ damage * pop + (1 | pop:fam)
##
## REML criterion at convergence: 2237.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2869 -0.6175 -0.0971 0.6606 3.9162
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 4.988 2.233
## Residual 33.217 5.763
## Number of obs: 351, groups: pop:fam, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 30.6963 1.0279 29.863
## damage1 1.9525 1.1440 1.707
## popteotihuacan 2.3377 1.3622 1.716
## popvaldeflores -0.1978 1.5611 -0.127
## popzubia 2.8334 1.7347 1.633
## damage1:popteotihuacan -1.7627 1.5413 -1.144
## damage1:popvaldeflores -1.8489 1.7887 -1.034
## damage1:popzubia -3.3676 1.9921 -1.691
##
## Correlation of Fixed Effects:
## (Intr) damag1 pptthc ppvldf popzub dmg1:ppt dmg1:ppv
## damage1 -0.596
## popteotihcn -0.755 0.450
## popvaldflrs -0.658 0.392 0.497
## popzubia -0.593 0.353 0.447 0.390
## dmg1:pptthc 0.442 -0.742 -0.588 -0.291 -0.262
## dmg1:ppvldf 0.381 -0.640 -0.288 -0.564 -0.226 0.475
## damag1:ppzb 0.342 -0.574 -0.258 -0.225 -0.564 0.426 0.367
##Modelo nulo
fit2<-lm(armf~1)
##Significancia modelo completo
anova(fit1, fit2)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## fit2: armf ~ 1
## fit1: armf ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 2280.5 2288.2 -1138.2 2276.5
## fit1 10 2273.3 2311.9 -1126.6 2253.3 23.208 8 0.003107 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia interacción
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## armf ~ damage * pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 2273.3
## damage:pop 3 2270.5 3.1838 0.3641
##Significancia del daño
fit3<-lmer(armf~pop+damage:pop+(1|pop:fam))
drop1(fit3, test = "Chisq")
## Single term deletions
##
## Model:
## armf ~ pop + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 2273.3
## pop:damage 4 2269.0 3.7333 0.4433
##Significancia de la población
fit4<-lmer(armf~damage+damage:pop+(1|pop:fam))
drop1(fit4, test = "Chisq")
## Single term deletions
##
## Model:
## armf ~ damage + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 2273.3
## damage:pop 6 2269.7 8.3822 0.2114
##Significancia efecto aleatorio fam
modsinfam<-lm(armf~damage*pop)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: armf ~ damage * pop
## fit1: armf ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 9 2281.6 2316.3 -1131.8 2263.6
## fit1 10 2273.3 2311.9 -1126.6 2253.3 10.312 1 0.001322 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Bootstrap paramétrico para reconfirmar significancia efecto del daño*fam (ho: la varianza es distinta de cero)
fit1f<-lmer(armf~damage*pop+(1|pop:fam), REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | pop:fam = 0
## from: lmer(formula = armf ~ damage * pop + (1 | pop:fam), REML = F)
## 1000 samples were taken Mon Oct 29 16:33:55 2018
## 3 samples had warnings, 3 in alternate model 0 in null model
## 3 unused samples. 0 <= P(abs((Intercept) | pop:fam) > |4.247412|) <= 0.003
#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$armffit2<-fit1fit
p<-ggplot(enemy,aes(damage,armffit2, 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)
##No.Hojas: Modelo general
fit1<-lmer(sqrt(leaves.2+0.5)~damage*pop+(1|pop:fam))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(leaves.2 + 0.5) ~ damage * pop + (1 | pop:fam)
##
## REML criterion at convergence: 277.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9760 -0.5390 -0.0720 0.4586 3.0523
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 0.04766 0.2183
## Residual 0.09908 0.3148
## Number of obs: 351, groups: pop:fam, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.76710 0.07421 37.289
## damage1 0.03957 0.06264 0.632
## popteotihuacan 0.17489 0.09814 1.782
## popvaldeflores -0.37544 0.11366 -3.303
## popzubia -0.35302 0.12660 -2.788
## damage1:popteotihuacan -0.06161 0.08441 -0.730
## damage1:popvaldeflores 0.15985 0.09787 1.633
## damage1:popzubia -0.02089 0.10920 -0.191
##
## Correlation of Fixed Effects:
## (Intr) damag1 pptthc ppvldf popzub dmg1:ppt dmg1:ppv
## damage1 -0.453
## popteotihcn -0.756 0.342
## popvaldflrs -0.653 0.295 0.494
## popzubia -0.586 0.265 0.443 0.383
## dmg1:pptthc 0.336 -0.742 -0.448 -0.219 -0.197
## dmg1:ppvldf 0.290 -0.640 -0.219 -0.424 -0.170 0.475
## damag1:ppzb 0.260 -0.574 -0.196 -0.169 -0.424 0.426 0.367
##Modelo nulo
fit2<-lm(sqrt(leaves.2+0.5)~1)
##Significancia modelo completo
anova(fit1, fit2)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## fit2: sqrt(leaves.2 + 0.5) ~ 1
## fit1: sqrt(leaves.2 + 0.5) ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 397.55 405.27 -196.77 393.55
## fit1 10 269.82 308.43 -124.91 249.82 143.72 8 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia interacción
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## sqrt(leaves.2 + 0.5) ~ damage * pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 269.82
## damage:pop 3 269.55 5.7226 0.1259
##Significancia del daño
fit3<-lmer(sqrt(leaves.2+0.5)~pop+damage:pop+(1|pop:fam))
drop1(fit3, test = "Chisq")
## Single term deletions
##
## Model:
## sqrt(leaves.2 + 0.5) ~ pop + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 269.82
## pop:damage 4 269.44 7.6117 0.1069
##Significancia de la población
fit4<-lmer(sqrt(leaves.2+0.5)~damage+damage:pop+(1|pop:fam))
drop1(fit4, test = "Chisq")
## Single term deletions
##
## Model:
## sqrt(leaves.2 + 0.5) ~ damage + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> 269.82
## damage:pop 6 290.73 32.908 1.092e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia efecto aleatorio fam
modsinfam<-lm(sqrt(leaves.2+0.5)~damage*pop)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: sqrt(leaves.2 + 0.5) ~ damage * pop
## fit1: sqrt(leaves.2 + 0.5) ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 9 318.07 352.82 -150.04 300.07
## fit1 10 269.82 308.43 -124.91 249.82 50.248 1 1.355e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Bootstrap paramétrico para reconfirmar significancia efecto del daño*fam (ho: la varianza es distinta de cero)
library(pbnm)
fit1f<-lmer(sqrt(leaves.2+0.5)~damage*pop+(1|pop:fam), REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | pop:fam = 0
## from: lmer(formula = sqrt(leaves.2 + 0.5) ~ damage * pop + (1 | pop:fam), REML = F)
## 1000 samples were taken Mon Oct 29 16:34:13 2018
## P(abs((Intercept) | pop:fam) > |0.04263713|) = 0
#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$leaves.2fit2<-fit1fit
library(ggplot2)
p<-ggplot(enemy,aes(damage,leaves.2fit2, 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)
##Masa total: Modelo general
fit1<-lmer(ln_wt.total~damage*pop+(1|pop:fam))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: ln_wt.total ~ damage * pop + (1 | pop:fam)
##
## REML criterion at convergence: -225.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2157 -0.5838 0.0090 0.4619 3.5420
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 0.01302 0.1141
## Residual 0.02241 0.1497
## Number of obs: 351, groups: pop:fam, 50
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.954812 0.037494 25.466
## damage1 -0.076758 0.029801 -2.576
## popteotihuacan 0.083460 0.049569 1.684
## popvaldeflores -0.101522 0.057504 -1.765
## popzubia -0.102834 0.064075 -1.605
## damage1:popteotihuacan 0.023691 0.040160 0.590
## damage1:popvaldeflores -0.009433 0.046560 -0.203
## damage1:popzubia -0.049688 0.051965 -0.956
##
## Correlation of Fixed Effects:
## (Intr) damag1 pptthc ppvldf popzub dmg1:ppt dmg1:ppv
## damage1 -0.426
## popteotihcn -0.756 0.322
## popvaldflrs -0.652 0.278 0.493
## popzubia -0.585 0.249 0.443 0.382
## dmg1:pptthc 0.316 -0.742 -0.422 -0.206 -0.185
## dmg1:ppvldf 0.273 -0.640 -0.206 -0.399 -0.160 0.475
## damag1:ppzb 0.244 -0.573 -0.185 -0.159 -0.399 0.426 0.367
##Modelo nulo
fit2<-lm(ln_wt.total~1)
##Significancia modelo completo
anova(fit1, fit2)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## fit2: ln_wt.total ~ 1
## fit1: ln_wt.total ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 -104.07 -96.352 54.037 -108.07
## fit1 10 -244.45 -205.846 132.227 -264.45 156.38 8 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia interacción
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## ln_wt.total ~ damage * pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> -244.45
## damage:pop 3 -248.24 2.2133 0.5293
##Significancia del daño
fit3<-lmer(ln_wt.total~pop+damage:pop+(1|pop:fam))
drop1(fit3, test = "Chisq")
## Single term deletions
##
## Model:
## ln_wt.total ~ pop + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> -244.45
## pop:damage 4 -228.04 24.413 6.601e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia de la población
fit4<-lmer(ln_wt.total~damage+damage:pop+(1|pop:fam))
drop1(fit4, test = "Chisq")
## Single term deletions
##
## Model:
## ln_wt.total ~ damage + damage:pop + (1 | pop:fam)
## Df AIC LRT Pr(Chi)
## <none> -244.45
## damage:pop 6 -233.05 23.402 0.0006724 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia efecto aleatorio fam
modsinfam<-lm(ln_wt.total~damage*pop)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: ln_wt.total ~ damage * pop
## fit1: ln_wt.total ~ damage * pop + (1 | pop:fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 9 -176.51 -141.77 97.257 -194.51
## fit1 10 -244.45 -205.85 132.227 -264.45 69.941 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Bootstrap paramétrico para reconfirmar significancia efecto del daño*fam (ho: la varianza es distinta de cero)
library(pbnm)
fit1f<-lmer(ln_wt.total~damage*pop+(1|pop:fam), REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | pop:fam = 0
## from: lmer(formula = ln_wt.total ~ damage * pop + (1 | pop:fam), REML = F)
## 1000 samples were taken Mon Oct 29 16:34:31 2018
## P(abs((Intercept) | pop:fam) > |0.01174032|) = 0
#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$ln_wt.totalfit2<-fit1fit
library(ggplot2)
p<-ggplot(enemy,aes(damage,ln_wt.totalfit2, 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)
##Guardando ajustados a la base
write.csv(enemy,file = "todos_2fit.csv")