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("valdeflores_fit.csv", sep = ",", header = T)
attach(enemy)
damage<-as.factor(damage)
fam<-as.factor(fam)
##ESFUERZO REPRODUCTIVO: Modelo general
fit1<-lmer(arep~damage+(1+damage|fam)+ln_wt.total)
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: arep ~ damage + (1 + damage | fam) + ln_wt.total
##
## REML criterion at convergence: 379.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4423 -0.5448 0.3117 0.6419 1.5483
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 5.428 2.330
## damage1 2.748 1.658 -1.00
## Residual 12.172 3.489
## Number of obs: 71, groups: fam, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 49.260 2.932 16.798
## damage1 -1.460 1.024 -1.426
## ln_wt.total -16.436 3.262 -5.039
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.499
## ln_wt.total -0.948 0.276
##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 + (1 + damage | fam) + ln_wt.total
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 421.71 426.24 -208.86 417.71
## fit1 7 400.12 415.96 -193.06 386.12 31.592 5 7.157e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia effectos fijos (daño y covariable)
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## arep ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC LRT Pr(Chi)
## <none> 400.12
## damage 1 400.18 2.0622 0.151
## ln_wt.total 1 419.45 21.3304 3.866e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia efecto aleatorio fam
modsinfam<-lm(arep~damage+ln_wt.total)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: arep ~ damage + ln_wt.total
## fit1: arep ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 4 398.40 407.45 -195.20 390.40
## fit1 7 400.12 415.96 -193.06 386.12 4.2792 3 0.2329
##Bootstrap paramétrico para reconfirmar significancia efecto del daño*fam (ho: la varianza es distinta de cero)
library(pbnm)
fit1f<-lmer(arep~damage+(1+damage|fam)+ln_wt.total, REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | fam damage | fam (Intercept)#damage | fam = 0
## from: lmer(formula = arep ~ damage + (1 + damage | fam) + ln_wt.total, REML = F)
## 1000 samples were taken Fri Oct 26 18:52:45 2018
## 133 samples had warnings, 133 in alternate model 0 in null model
## 133 unused samples. 0.047 <= P(LRT stat > |4.279161|) <= 0.18
#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$arepfit<-fit1fit
library(ggplot2)
p<-ggplot(enemy,aes(damage,arepfit, group=fam))
p + geom_smooth(method = "lm", se=F)
##LMF: Modelo general
fit1<-lmer(almf~damage+(1+damage|fam)+ln_wt.total)
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: almf ~ damage + (1 + damage | fam) + ln_wt.total
##
## REML criterion at convergence: 394.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4083 -0.4852 -0.1576 0.3800 6.0322
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 2.3164 1.5220
## damage1 0.8965 0.9468 -1.00
## Residual 16.3939 4.0489
## Number of obs: 71, groups: fam, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 17.0704 3.1895 5.352
## damage1 -0.1952 1.0586 -0.184
## ln_wt.total 0.7538 3.6135 0.209
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.460
## ln_wt.total -0.967 0.299
##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 + (1 + damage | fam) + ln_wt.total
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 406.50 411.03 -201.25 402.50
## fit1 7 415.66 431.50 -200.83 401.66 0.835 5 0.9747
##Significancia effectos fijos (daño y covariable)
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## almf ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC LRT Pr(Chi)
## <none> 415.66
## damage 1 413.70 0.035290 0.8510
## ln_wt.total 1 413.71 0.045364 0.8313
##Significancia efecto aleatorio fam
modsinfam<-lm(almf~damage+ln_wt.total)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: almf ~ damage + ln_wt.total
## fit1: almf ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 4 410.38 419.43 -201.19 402.38
## fit1 7 415.66 431.50 -200.83 401.66 0.7167 3 0.8693
##Bootstrap paramétrico para reconfirmar significancia efecto del daño*fam (ho: la varianza es distinta de cero)
fit1f<-lmer(almf~damage+(1+damage|fam)+ln_wt.total, REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | fam damage | fam (Intercept)#damage | fam = 0
## from: lmer(formula = almf ~ damage + (1 + damage | fam) + ln_wt.total, REML = F)
## 1000 samples were taken Fri Oct 26 18:53:04 2018
## 133 samples had warnings, 133 in alternate model 0 in null model
## 133 unused samples. 0.33 <= P(LRT stat > |0.7167152|) <= 0.463
#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$almffit<-fit1fit
p<-ggplot(enemy,aes(damage,almffit, group=fam))
p + geom_smooth(method = "lm", se=F)
##SMF: Modelo general
fit1<-lmer(asmf~damage+(1+damage|fam)+ln_wt.total)
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: asmf ~ damage + (1 + damage | fam) + ln_wt.total
##
## REML criterion at convergence: 365.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5490 -0.6300 -0.0838 0.5781 3.2452
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 3.3495 1.8302
## damage1 0.4239 0.6511 -1.00
## Residual 9.9349 3.1520
## Number of obs: 71, groups: fam, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 34.8604 2.6570 13.120
## damage1 0.2804 0.8214 0.341
## ln_wt.total -1.6971 2.9794 -0.570
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.476
## ln_wt.total -0.957 0.315
##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 + (1 + damage | fam) + ln_wt.total
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 379.44 383.97 -187.72 375.44
## fit1 7 385.43 401.27 -185.72 371.43 4.0075 5 0.5483
##Significancia effectos fijos (daño y covariable)
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## asmf ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC LRT Pr(Chi)
## <none> 385.43
## damage 1 383.56 0.12488 0.7238
## ln_wt.total 1 383.67 0.23780 0.6258
##Significancia efecto aleatorio fam
modsinfam<-lm(asmf~damage+ln_wt.total)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: asmf ~ damage + ln_wt.total
## fit1: asmf ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 4 383.25 392.30 -187.62 375.25
## fit1 7 385.43 401.27 -185.72 371.43 3.8125 3 0.2824
##Bootstrap paramétrico para reconfirmar significancia efecto del daño*fam (ho: la varianza es distinta de cero)
fit1f<-lmer(asmf~damage+(1+damage|fam)+ln_wt.total, REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | fam damage | fam (Intercept)#damage | fam = 0
## from: lmer(formula = asmf ~ damage + (1 + damage | fam) + ln_wt.total, REML = F)
## 1000 samples were taken Fri Oct 26 18:53:23 2018
## 136 samples had warnings, 136 in alternate model 0 in null model
## 136 unused samples. 0.056 <= P(LRT stat > |3.812548|) <= 0.192
#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$asmffit<-fit1fit
p<-ggplot(enemy,aes(damage,asmffit, group=fam))
p + geom_smooth(method = "lm", se=F)
##RMF: Modelo general
fit1<-lmer(armf~damage+(1+damage|fam)+ln_wt.total)
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: armf ~ damage + (1 + damage | fam) + ln_wt.total
##
## REML criterion at convergence: 425.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.47131 -0.67794 0.01606 0.62843 2.70606
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 3.3099 1.8193
## damage1 0.1962 0.4429 -1.00
## Residual 25.6506 5.0646
## Number of obs: 71, groups: fam, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 16.893 4.040 4.181
## damage1 1.525 1.280 1.191
## ln_wt.total 15.934 4.585 3.475
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.450
## ln_wt.total -0.969 0.313
##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 + (1 + damage | fam) + ln_wt.total
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 451.67 456.20 -223.84 447.67
## fit1 7 448.22 464.06 -217.11 434.22 13.448 5 0.01953 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia effectos fijos (daño y covariable)
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## armf ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC LRT Pr(Chi)
## <none> 448.22
## damage 1 447.68 1.4593 0.2270387
## ln_wt.total 1 457.72 11.4909 0.0006994 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia efecto aleatorio fam
modsinfam<-lm(armf~damage+ln_wt.total)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: armf ~ damage + ln_wt.total
## fit1: armf ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 4 443.31 452.36 -217.65 435.31
## fit1 7 448.22 464.06 -217.11 434.22 1.0838 3 0.781
##Bootstrap paramétrico para reconfirmar significancia efecto del daño*fam (ho: la varianza es distinta de cero)
fit1f<-lmer(armf~damage+(1+damage|fam)+ln_wt.total, REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | fam damage | fam (Intercept)#damage | fam = 0
## from: lmer(formula = armf ~ damage + (1 + damage | fam) + ln_wt.total, REML = F)
## 1000 samples were taken Fri Oct 26 18:53:41 2018
## 136 samples had warnings, 136 in alternate model 0 in null model
## 136 unused samples. 0.254 <= P(LRT stat > |1.083752|) <= 0.39
#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$armffit<-fit1fit
p<-ggplot(enemy,aes(damage,armffit, group=fam))
p + geom_smooth(method = "lm", se=F)
##No.Hojas: Modelo general
fit1<-lmer(sqrt(leaves.2+0.5)~damage+(1+damage|fam)+ln_wt.total)
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(leaves.2 + 0.5) ~ damage + (1 + damage | fam) + ln_wt.total
##
## REML criterion at convergence: 27.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3803 -0.5885 -0.0716 0.6084 2.5728
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 0.01390 0.1179
## damage1 0.02476 0.1573 -0.16
## Residual 0.06513 0.2552
## Number of obs: 71, groups: fam, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.84342 0.21712 8.490
## damage1 0.25349 0.08177 3.100
## ln_wt.total 0.64503 0.24593 2.623
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.363
## ln_wt.total -0.966 0.257
##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 + (1 + damage | fam) + ln_wt.total
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 44.416 48.942 -20.2082 40.416
## fit1 7 33.213 49.052 -9.6065 19.213 21.203 5 0.0007414 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia effectos fijos (daño y covariable)
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## sqrt(leaves.2 + 0.5) ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC LRT Pr(Chi)
## <none> 33.213
## damage 1 38.683 7.4703 0.006273 **
## ln_wt.total 1 38.027 6.8142 0.009043 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia efecto aleatorio fam
modsinfam<-lm(sqrt(leaves.2+0.5)~damage+ln_wt.total)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: sqrt(leaves.2 + 0.5) ~ damage + ln_wt.total
## fit1: sqrt(leaves.2 + 0.5) ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 4 32.746 41.797 -12.3730 24.746
## fit1 7 33.213 49.052 -9.6065 19.213 5.5331 3 0.1367
##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+(1+damage|fam)+ln_wt.total, REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | fam damage | fam (Intercept)#damage | fam = 0
## from: lmer(formula = sqrt(leaves.2 + 0.5) ~ damage + (1 + damage | fam) + ln_wt.total, REML = F)
## 1000 samples were taken Fri Oct 26 18:54:00 2018
## 134 samples had warnings, 134 in alternate model 0 in null model
## 134 unused samples. 0.022 <= P(LRT stat > |5.533075|) <= 0.156
#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.2fit<-fit1fit
library(ggplot2)
p<-ggplot(enemy,aes(damage,leaves.2fit, group=fam))
p + geom_smooth(method = "lm", se=F)
##Masa total: Modelo general
fit1<-lmer(ln_wt.total~damage+(1+damage|fam))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: ln_wt.total ~ damage + (1 + damage | fam)
##
## REML criterion at convergence: -72.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.37633 -0.52729 0.06114 0.53952 2.79684
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 0.0064305 0.08019
## damage1 0.0006402 0.02530 -0.61
## Residual 0.0155154 0.12456
## Number of obs: 71, groups: fam, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.85310 0.03253 26.227
## damage1 -0.08640 0.03082 -2.803
##
## Correlation of Fixed Effects:
## (Intr)
## damage1 -0.533
##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 + (1 + damage | fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit2 2 -65.147 -60.622 34.573 -69.147
## fit1 6 -71.426 -57.850 41.713 -83.426 14.279 4 0.006456 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia effectos fijos (daño y covariable)
drop1(fit1, test = "Chisq")
## Single term deletions
##
## Model:
## ln_wt.total ~ damage + (1 + damage | fam)
## Df AIC LRT Pr(Chi)
## <none> -71.426
## damage 1 -67.229 6.1971 0.0128 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Significancia efecto aleatorio fam
modsinfam<-lm(ln_wt.total~damage)
anova(fit1, modsinfam)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## modsinfam: ln_wt.total ~ damage
## fit1: ln_wt.total ~ damage + (1 + damage | fam)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modsinfam 3 -69.695 -62.907 37.847 -75.695
## fit1 6 -71.426 -57.850 41.713 -83.426 7.7315 3 0.0519 .
## ---
## 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+(1+damage|fam), REML = F)
bootfam<-pbnm(fit1f,modsinfam,nsim = 1000, tasks = 10, cores=2, seed=3400835)
summary(bootfam)
## Parametric bootstrap testing: (Intercept) | fam damage | fam (Intercept)#damage | fam = 0
## from: lmer(formula = ln_wt.total ~ damage + (1 + damage | fam), REML = F)
## 1000 samples were taken Fri Oct 26 18:54:19 2018
## 138 samples had warnings, 138 in alternate model 0 in null model
## 138 unused samples. 0.01 <= P(LRT stat > |7.731461|) <= 0.148
#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.totalfit<-fit1fit
library(ggplot2)
p<-ggplot(enemy,aes(damage,ln_wt.totalfit, group=fam))
p + geom_smooth(method = "lm", se=F)
##Guardando ajustados a la base
write.csv(enemy, file = "valdeflores_fit2.csv")