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("teo_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: 742
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5451 -0.4590 0.1332 0.5845 3.2713
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 1.888 1.374
## damage1 1.751 1.323 1.00
## Residual 19.693 4.438
## Number of obs: 126, groups: fam, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 43.4336 2.3245 18.685
## damage1 -0.5099 0.8592 -0.593
## ln_wt.total -10.9046 2.1491 -5.074
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.237
## ln_wt.total -0.961 0.132
##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 785.22 790.89 -390.61 781.22
## fit1 7 761.62 781.48 -373.81 747.62 33.598 5 2.863e-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> 761.62
## damage 1 759.98 0.3538 0.5519
## ln_wt.total 1 782.63 23.0063 1.615e-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 763.42 774.76 -377.71 755.42
## fit1 7 761.62 781.48 -373.81 747.62 7.7952 3 0.05044 .
## ---
## 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:00:21 2018
## 122 samples had warnings, 122 in alternate model 0 in null model
## 122 unused samples. 0.015 <= P(LRT stat > |7.795197|) <= 0.137
#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: 676.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7382 -0.3803 -0.0222 0.3647 6.5513
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 3.398 1.843
## damage1 6.991 2.644 -1.00
## Residual 11.863 3.444
## Number of obs: 126, groups: fam, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 24.8091 1.6261 15.256
## damage1 -0.6373 0.8707 -0.732
## ln_wt.total -6.5947 1.4463 -4.560
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.399
## ln_wt.total -0.927 0.087
##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 705.06 710.73 -350.53 701.06
## fit1 7 694.39 714.24 -340.19 680.39 20.67 5 0.0009351 ***
## ---
## 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> 694.39
## damage 1 692.94 0.5559 0.4559
## ln_wt.total 1 710.57 18.1803 2.009e-05 ***
## ---
## 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 692.89 704.23 -342.44 684.89
## fit1 7 694.39 714.24 -340.19 680.39 4.5001 3 0.2123
##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:00:41 2018
## 122 samples had warnings, 122 in alternate model 0 in null model
## 122 unused samples. 0.05 <= P(LRT stat > |4.500092|) <= 0.172
#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: 731.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9431 -0.5232 -0.0837 0.4527 3.5931
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 11.6422 3.4121
## damage1 0.4876 0.6983 -1.00
## Residual 16.4904 4.0608
## Number of obs: 126, groups: fam, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 39.8028 2.4139 16.489
## damage1 -0.2411 0.7567 -0.319
## ln_wt.total -5.7458 2.1408 -2.684
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.355
## ln_wt.total -0.921 0.150
##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 763.59 769.26 -379.79 759.59
## fit1 7 751.31 771.16 -368.65 737.31 22.281 5 0.000463 ***
## ---
## 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> 751.31
## damage 1 749.41 0.0978 0.75455
## ln_wt.total 1 755.59 6.2832 0.01219 *
## ---
## 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 766.11 777.45 -379.05 758.11
## fit1 7 751.31 771.16 -368.65 737.31 20.798 3 0.000116 ***
## ---
## 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:01:01 2018
## 122 samples had warnings, 122 in alternate model 0 in null model
## 122 unused samples. 0 <= P(LRT stat > |20.79771|) <= 0.122
#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: 768.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0988 -0.5093 0.0497 0.5392 2.9458
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 5.592 2.365
## damage1 2.219 1.490 -1.00
## Residual 25.607 5.060
## Number of obs: 126, groups: fam, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 15.7928 2.5348 6.230
## damage1 1.0989 0.9775 1.124
## ln_wt.total 16.5671 2.2975 7.211
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.360
## ln_wt.total -0.943 0.123
##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 821.99 827.66 -409.0 817.99
## fit1 7 788.41 808.26 -387.2 774.41 43.581 5 2.818e-08 ***
## ---
## 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> 788.41
## damage 1 787.62 1.215 0.2703
## ln_wt.total 1 828.44 42.035 8.967e-11 ***
## ---
## 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 784.68 796.02 -388.34 776.68
## fit1 7 788.41 808.26 -387.20 774.41 2.2675 3 0.5188
##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:01:21 2018
## 124 samples had warnings, 124 in alternate model 0 in null model
## 124 unused samples. 0.161 <= P(LRT stat > |2.267515|) <= 0.285
#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: 121
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9538 -0.4604 0.0119 0.3991 2.2508
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 0.03074 0.1753
## damage1 0.02253 0.1501 0.42
## Residual 0.11653 0.3414
## Number of obs: 126, groups: fam, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.09200 0.19282 16.036
## damage1 -0.02696 0.07107 -0.379
## ln_wt.total -0.14721 0.17659 -0.834
##
## Correlation of Fixed Effects:
## (Intr) damag1
## damage1 -0.221
## ln_wt.total -0.952 0.133
##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 131.11 136.79 -63.557 127.11
## fit1 7 125.89 145.74 -55.944 111.89 15.226 5 0.009441 **
## ---
## 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> 125.89
## damage 1 124.03 0.14252 0.7058
## ln_wt.total 1 124.50 0.60720 0.4358
##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 134.03 145.38 -63.016 126.03
## fit1 7 125.89 145.74 -55.944 111.89 14.143 3 0.002717 **
## ---
## 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:01:40 2018
## 120 samples had warnings, 120 in alternate model 0 in null model
## 120 unused samples. 0 <= P(LRT stat > |14.14263|) <= 0.12
#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: -52.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.89732 -0.53107 0.00964 0.44876 2.95297
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## fam (Intercept) 2.050e-02 0.143178
## damage1 5.448e-05 0.007381 1.00
## Residual 2.758e-02 0.166060
## Number of obs: 126, groups: fam, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.03804 0.03911 26.541
## damage1 -0.05309 0.02992 -1.774
##
## Correlation of Fixed Effects:
## (Intr)
## damage1 -0.335
##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 -19.296 -13.624 11.648 -23.296
## fit1 6 -50.517 -33.500 31.259 -62.517 39.221 4 6.271e-08 ***
## ---
## 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> -50.517
## damage 1 -49.387 3.1303 0.07685 .
## ---
## 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 -19.059 -10.551 12.530 -25.059
## fit1 6 -50.517 -33.500 31.259 -62.517 37.458 3 3.681e-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)
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:02:01 2018
## 104 samples had warnings, 104 in alternate model 0 in null model
## 104 unused samples. 0 <= P(LRT stat > |37.45804|) <= 0.104
#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 = "teo_fit2.csv")