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("invas_fam_out_2.csv", sep = ",", header = T)
attach(enemy)
damage<-as.factor(damage)
fam<-as.factor(fam)
##Modelo sin interacciónes
fit1<-lmer(ln.wt.leaves.2~damage*origen+(1|pop))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: ln.wt.leaves.2 ~ damage * origen + (1 | pop)
##
## REML criterion at convergence: -756.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7002 -0.3985 -0.0977 0.2405 11.1439
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop (Intercept) 3.694e-17 6.078e-09
## Residual 6.859e-03 8.282e-02
## Number of obs: 365, groups: pop, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.11874 0.01019 11.647
## damage1 -0.01470 0.01491 -0.986
## origenMEX 0.04354 0.01275 3.415
## damage1:origenMEX -0.01649 0.01833 -0.899
##
## Correlation of Fixed Effects:
## (Intr) damag1 orgMEX
## damage1 -0.684
## origenMEX -0.800 0.547
## dmg1:rgnMEX 0.556 -0.813 -0.695
##Probando efectos fijos (damage) con la correción Kenward-Roger para datos desbalanceados
library(pbkrtest)
modtodos<-lmer(ln.wt.leaves.2~damage+origen+(1|pop)+damage:origen, REML = F)
modsindamage<-lmer(ln.wt.leaves.2~origen+(1|pop), REML = F)
p.damage<-KRmodcomp(modtodos,modsindamage)
p.damage
## F-test with Kenward-Roger approximation; computing time: 0.14 sec.
## large : ln.wt.leaves.2 ~ damage + origen + (1 | pop) + damage:origen
## small : ln.wt.leaves.2 ~ origen + (1 | pop)
## stat ndf ddf F.scaling p.value
## Ftest 4.7535 2.0000 359.2207 1 0.009171 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modsinorigen<-lmer(ln.wt.leaves.2~damage+(1|pop), REML = F)
p.origen<-KRmodcomp(modtodos,modsinorigen)
p.origen
## F-test with Kenward-Roger approximation; computing time: 0.08 sec.
## large : ln.wt.leaves.2 ~ damage + origen + (1 | pop) + damage:origen
## small : ln.wt.leaves.2 ~ damage + (1 | pop)
## stat ndf ddf F.scaling p.value
## Ftest 7.0259 2.0000 6.1867 0.90455 0.02557 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modsininter<-lmer(ln.wt.leaves.2~damage+origen+(1|pop), REML=F)
p.inter<-KRmodcomp(modtodos,modsininter)
p.inter
## F-test with Kenward-Roger approximation; computing time: 0.08 sec.
## large : ln.wt.leaves.2 ~ damage + origen + (1 | pop) + damage:origen
## small : ln.wt.leaves.2 ~ damage + origen + (1 | pop)
## stat ndf ddf F.scaling p.value
## Ftest 0.8089 1.0000 359.2169 1 0.369
##Probando efectos aleatorios
library(lmerTest)
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
ranova(fit1)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## ln.wt.leaves.2 ~ damage + origen + (1 | pop) + damage:origen
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 6 378.14 -744.27
## (1 | pop) 5 378.14 -746.27 1.1369e-13 1 1
detach("package:lmerTest", unload=TRUE)
#Finalmente el diagnostico de residuales del modelo...
plot(residuals(fit1)~predict(fit1,type="link"), xlab=expression(hat(eta)), ylab="Deviance residuals")
#...y los intervalos de confianza para las familias y poblaciones
library(lattice)
ranef(fit1)
## $pop
## (Intercept)
## morelos -4.109625e-15
## teotihuacan 4.109625e-15
## valdeflores 3.616223e-16
## zubia -3.616223e-16
dput(cl<-ranef(fit1,condVar=T))
## structure(list(pop = structure(list(`(Intercept)` = c(-4.10962512548812e-15,
## 4.10962512548816e-15, 3.61622319201582e-16, -3.61622319201281e-16
## )), class = "data.frame", row.names = c("morelos", "teotihuacan",
## "valdeflores", "zubia"), postVar = structure(c(3.69368739979627e-17,
## 3.69368739979593e-17, 3.69368739979704e-17, 3.69368739979748e-17
## ), .Dim = c(1L, 1L, 4L)))), class = "ranef.mer")
dotplot(cl)
## $pop
