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("morelos_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: 605
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.62406 -0.58909 -0.01562 0.49509 2.42272
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 4.74892 2.1792
## damage1 0.04378 0.2092 1.00
## Residual 19.80942 4.4508
## Number of obs: 103, groups: fam, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 42.8988 3.1005 13.836
## damage1 -1.1608 0.9181 -1.264
## ln_wt.total -9.9406 3.1149 -3.191
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.386
## ln_wt.total -0.960 0.261
##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 638.75 644.02 -317.37 634.75
## fit1 7 625.92 644.36 -305.96 611.92 22.833 5 0.0003633 ***
## ---
## 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> 625.92
## damage 1 625.56 1.6419 0.200059
## ln_wt.total 1 633.69 9.7755 0.001769 **
## ---
## 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 627.11 637.65 -309.56 619.11
## fit1 7 625.92 644.36 -305.96 611.92 7.1986 3 0.06583 .
## ---
## 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+(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:41:30 2018
## 130 samples had warnings, 130 in alternate model 0 in null model
## 130 unused samples. 0.015 <= P(LRT stat > |7.198625|) <= 0.145
#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: 597.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6466 -0.4297 -0.0246 0.2700 7.1824
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 3.724 1.930
## damage1 2.046 1.431 -1.00
## Residual 19.774 4.447
## Number of obs: 103, groups: fam, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 22.8954 2.8740 7.967
## damage1 -1.4532 0.9873 -1.472
## ln_wt.total -4.2634 2.8791 -1.481
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.432
## ln_wt.total -0.958 0.224
##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 614.99 620.26 -305.50 610.99
## fit1 7 617.95 636.40 -301.98 603.95 7.0411 5 0.2176
##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> 617.95
## damage 1 618.15 2.1930 0.1386
## ln_wt.total 1 618.01 2.0599 0.1512
##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 613.16 623.7 -302.58 605.16
## fit1 7 617.95 636.4 -301.98 603.95 1.2119 3 0.7502
##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:41:50 2018
## 132 samples had warnings, 132 in alternate model 0 in null model
## 132 unused samples. 0.267 <= P(LRT stat > |1.211905|) <= 0.399
#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: 576.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2383 -0.5295 0.1100 0.6035 1.8524
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 7.8788 2.8069
## damage1 0.1722 0.4149 -1.00
## Residual 14.0544 3.7489
## Number of obs: 103, groups: fam, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 39.9143 2.7584 14.470
## damage1 -0.8218 0.7828 -1.050
## ln_wt.total -6.1094 2.7204 -2.246
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.428
## ln_wt.total -0.942 0.266
##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 600.89 606.16 -298.45 596.89
## fit1 7 596.50 614.94 -291.25 582.50 14.398 5 0.01327 *
## ---
## 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:
## asmf ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC LRT Pr(Chi)
## <none> 596.50
## damage 1 595.56 1.0666 0.30172
## ln_wt.total 1 598.59 4.0943 0.04303 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##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 604.6 615.13 -298.30 596.6
## fit1 7 596.5 614.94 -291.25 582.5 14.099 3 0.002774 **
## ---
## 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+(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:42:09 2018
## 133 samples had warnings, 133 in alternate model 0 in null model
## 133 unused samples. 0 <= P(LRT stat > |14.0987|) <= 0.133
#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: 640.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.15974 -0.61191 -0.07557 0.57473 2.63301
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 7.3305 2.7075
## damage1 0.1351 0.3676 -1.00
## Residual 28.6223 5.3500
## Number of obs: 103, groups: fam, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 18.179 3.702 4.910
## damage1 2.949 1.105 2.668
## ln_wt.total 13.104 3.712 3.530
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.413
## ln_wt.total -0.958 0.258
##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 675.21 680.48 -335.6 671.21
## fit1 7 662.19 680.64 -324.1 648.19 23.013 5 0.0003356 ***
## ---
## 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> 662.19
## damage 1 667.06 6.8613 0.008808 **
## ln_wt.total 1 672.09 11.8978 0.000562 ***
## ---
## 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 661.79 672.33 -326.89 653.79
## fit1 7 662.19 680.64 -324.10 648.19 5.5933 3 0.1332
##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:42:29 2018
## 132 samples had warnings, 132 in alternate model 0 in null model
## 132 unused samples. 0.03 <= P(LRT stat > |5.593293|) <= 0.162
#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: 97.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.13019 -0.75839 -0.08127 0.52843 2.79716
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 0.036778 0.19178
## damage1 0.003363 0.05799 1.00
## Residual 0.118364 0.34404
## Number of obs: 103, groups: fam, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.93399 0.24558 11.947
## damage1 0.02531 0.07265 0.348
## ln_wt.total -0.17300 0.24593 -0.703
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.345
## ln_wt.total -0.957 0.261
##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 108.51 113.78 -52.255 104.510
## fit1 7 103.25 121.69 -44.623 89.246 15.264 5 0.009291 **
## ---
## 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> 103.25
## damage 1 101.36 0.11287 0.7369
## ln_wt.total 1 101.78 0.53233 0.4656
##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 110.22 120.76 -51.112 102.223
## fit1 7 103.25 121.69 -44.623 89.246 12.977 3 0.004685 **
## ---
## 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+(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:42:48 2018
## 132 samples had warnings, 132 in alternate model 0 in null model
## 132 unused samples. 0.001 <= P(LRT stat > |12.97753|) <= 0.133
#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: -87.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7075 -0.5922 -0.0494 0.5146 3.2998
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 0.0115585 0.10751
## damage1 0.0004127 0.02032 -1.00
## Residual 0.0186965 0.13674
## Number of obs: 103, groups: fam, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.95477 0.03499 27.287
## damage1 -0.07651 0.02777 -2.755
##
## Correlation of Fixed Effects:
## (Intr)
## damage1 -0.571
##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 -71.983 -66.714 37.992 -75.983
## fit1 6 -85.992 -70.183 48.996 -97.992 22.009 4 0.0001996 ***
## ---
## 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> -85.992
## damage 1 -80.882 7.1092 0.007669 **
## ---
## 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 -75.697 -67.793 40.849 -81.697
## fit1 6 -85.992 -70.183 48.996 -97.992 16.294 3 0.0009867 ***
## ---
## 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:43:08 2018
## 125 samples had warnings, 125 in alternate model 0 in null model
## 125 unused samples. 0 <= P(LRT stat > |16.29447|) <= 0.125
#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 = "morelos_fit2.csv")