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("zubia_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: 269.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0781 -0.6029 0.1751 0.5317 1.9095
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 4.918 2.218
## damage1 4.558 2.135 -0.71
## Residual 11.679 3.417
## Number of obs: 51, groups: fam, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 49.488 2.889 17.133
## damage1 -1.939 1.318 -1.471
## ln_wt.total -17.954 3.146 -5.707
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.512
## ln_wt.total -0.929 0.290
##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 307.60 311.46 -151.80 303.60
## fit1 7 291.29 304.81 -138.64 277.29 26.309 5 7.773e-05 ***
## ---
## 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> 291.29
## damage 1 291.54 2.2558 0.1331
## ln_wt.total 1 313.74 24.4483 7.633e-07 ***
## ---
## 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 287.82 295.55 -139.91 279.82
## fit1 7 291.29 304.81 -138.64 277.29 2.5325 3 0.4695
##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:47:33 2018
## 158 samples had warnings, 158 in alternate model 0 in null model
## 158 unused samples. 0.104 <= P(LRT stat > |2.532453|) <= 0.262
#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: 223.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.70368 -0.47330 -0.08688 0.60854 2.16152
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 3.477 1.865
## damage1 3.084 1.756 -0.89
## Residual 4.307 2.075
## Number of obs: 51, groups: fam, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 24.3314 1.8610 13.074
## damage1 0.7104 0.9190 0.773
## ln_wt.total -8.1290 1.9591 -4.149
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.574
## ln_wt.total -0.900 0.261
##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 255.95 259.81 -125.98 251.95
## fit1 7 242.81 256.33 -114.40 228.81 23.146 5 0.0003166 ***
## ---
## 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:
## almf ~ damage + (1 + damage | fam) + ln_wt.total
## Df AIC LRT Pr(Chi)
## <none> 242.81
## damage 1 241.51 0.702 0.4021229
## ln_wt.total 1 253.72 12.911 0.0003267 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##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 242.54 250.27 -117.27 234.54
## fit1 7 242.81 256.33 -114.40 228.81 5.7378 3 0.1251
##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:47:52 2018
## 158 samples had warnings, 158 in alternate model 0 in null model
## 158 unused samples. 0.023 <= P(LRT stat > |5.737807|) <= 0.181
#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: 251.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.41341 -0.55778 0.03291 0.57370 3.07222
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 0.000e+00 0.000e+00
## damage1 4.335e-15 6.584e-08 NaN
## Residual 9.616e+00 3.101e+00
## Number of obs: 51, groups: fam, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 33.41107 2.15761 15.485
## damage1 -0.06607 0.91257 -0.072
## ln_wt.total -1.86162 2.45065 -0.760
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.471
## ln_wt.total -0.961 0.302
##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 261.71 265.57 -128.85 257.71
## fit1 7 271.07 284.60 -128.54 257.07 0.638 5 0.9862
##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> 271.07
## damage 1 269.08 0.00557 0.9405
## ln_wt.total 1 269.68 0.60947 0.4350
##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 265.07 272.8 -128.54 257.07
## fit1 7 271.07 284.6 -128.54 257.07 0 3 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:48:10 2018
## 157 samples had warnings, 157 in alternate model 0 in null model
## 157 unused samples. 0.565 <= P(LRT stat > |1.136868e-13|) <= 0.722
#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: 292.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.67199 -0.65267 0.03642 0.58694 2.03429
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 1.25381 1.1197
## damage1 0.09763 0.3125 -1.00
## Residual 21.50435 4.6373
## Number of obs: 51, groups: fam, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 14.425 3.382 4.265
## damage1 1.296 1.385 0.936
## ln_wt.total 22.488 3.813 5.898
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.490
## ln_wt.total -0.956 0.322
##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 333.05 336.92 -164.53 329.05
## fit1 7 314.02 327.54 -150.01 300.02 29.034 5 2.284e-05 ***
## ---
## 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> 314.02
## damage 1 312.91 0.8862 0.3465
## ln_wt.total 1 339.77 27.7500 1.381e-07 ***
## ---
## 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 308.05 315.77 -150.02 300.05
## fit1 7 314.02 327.54 -150.01 300.02 0.0248 3 0.999
##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:48:28 2018
## 155 samples had warnings, 155 in alternate model 0 in null model
## 155 unused samples. 0.494 <= P(LRT stat > |0.02476092|) <= 0.649
#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: -4.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.29727 -0.41229 0.04363 0.48046 2.10173
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 0.084819 0.29124
## damage1 0.003887 0.06235 -1.00
## Residual 0.033028 0.18174
## Number of obs: 51, groups: fam, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.04750 0.18612 11.001
## damage1 0.07317 0.06100 1.200
## ln_wt.total 0.43063 0.17077 2.522
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.624
## ln_wt.total -0.784 0.363
##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 28.5626 32.426 -12.2813 24.563
## fit1 7 1.2154 14.738 6.3923 -12.785 37.347 5 5.102e-07 ***
## ---
## 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> 1.2154
## damage 1 0.7117 1.4963 0.22125
## ln_wt.total 1 5.1907 5.9752 0.01451 *
## ---
## 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 29.8068 37.534 -10.9034 21.807
## fit1 7 1.2154 14.738 6.3923 -12.785 34.591 3 1.486e-07 ***
## ---
## 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:48:47 2018
## 157 samples had warnings, 157 in alternate model 0 in null model
## 157 unused samples. 0 <= P(LRT stat > |34.59133|) <= 0.157
#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: -28.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6688 -0.5410 -0.1442 0.3887 3.2781
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 0.01685 0.1298
## damage1 0.01414 0.1189 -0.86
## Residual 0.02301 0.1517
## Number of obs: 51, groups: fam, 7
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.85445 0.05725 14.926
## damage1 -0.12277 0.06221 -1.973
##
## Correlation of Fixed Effects:
## (Intr)
## damage1 -0.778
##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 -22.906 -19.043 13.453 -26.906
## fit1 6 -25.277 -13.686 18.638 -37.277 10.371 4 0.03463 *
## ---
## 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> -25.277
## damage 1 -23.880 3.397 0.06531 .
## ---
## 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 -25.786 -19.991 15.893 -31.786
## fit1 6 -25.277 -13.686 18.638 -37.277 5.4904 3 0.1392
##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:49:05 2018
## 166 samples had warnings, 166 in alternate model 0 in null model
## 166 unused samples. 0.027 <= P(LRT stat > |5.490397|) <= 0.193
#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 = "zubia_fit2.csv")