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(rel.fitness~damage*pop+(1|pop:fam))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: rel.fitness ~ damage * pop + (1 | pop:fam)
##
## REML criterion at convergence: 450.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4568 -0.4069 0.0587 0.4900 5.0123
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 0.0235 0.1533
## Residual 0.1675 0.4093
## Number of obs: 376, groups: pop:fam, 52
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.17166 0.06765 17.320
## damage1 -0.13922 0.07647 -1.821
## popteotihuacan -0.06628 0.09149 -0.724
## popvaldeflores -0.32151 0.10615 -3.029
## popzubia -0.26646 0.11910 -2.237
## damage1:popteotihuacan 0.10892 0.10440 1.043
## damage1:popvaldeflores 0.05262 0.12233 0.430
## damage1:popzubia -0.05422 0.13875 -0.391
##
## Correlation of Fixed Effects:
## (Intr) damag1 pptthc ppvldf popzub dmg1:ppt dmg1:ppv
## damage1 -0.584
## popteotihcn -0.739 0.432
## popvaldflrs -0.637 0.372 0.471
## popzubia -0.568 0.332 0.420 0.362
## dmg1:pptthc 0.428 -0.732 -0.589 -0.272 -0.243
## dmg1:ppvldf 0.365 -0.625 -0.270 -0.565 -0.207 0.458
## damag1:ppzb 0.322 -0.551 -0.238 -0.205 -0.564 0.404 0.344
##Probando efectos fijos (damage) con la correción Kenward-Roger para datos desbalanceados
library(pbkrtest)
modtodos<-lmer(rel.fitness~damage+pop+(1|pop:fam)+damage:pop, REML = F)
modsindamage<-lmer(rel.fitness~pop+(1|pop:fam), REML = F)
p.damage<-KRmodcomp(modtodos,modsindamage)
p.damage
## F-test with Kenward-Roger approximation; computing time: 0.12 sec.
## large : rel.fitness ~ damage + pop + (1 | pop:fam) + damage:pop
## small : rel.fitness ~ pop + (1 | pop:fam)
## stat ndf ddf F.scaling p.value
## Ftest 1.7757 4.0000 328.3464 1 0.1334
modsinpop<-lmer(rel.fitness~damage+(1|pop:fam), REML = F)
p.pop<-KRmodcomp(modtodos,modsinpop)
p.pop
## F-test with Kenward-Roger approximation; computing time: 0.07 sec.
## large : rel.fitness ~ damage + pop + (1 | pop:fam) + damage:pop
## small : rel.fitness ~ damage + (1 | pop:fam)
## stat ndf ddf F.scaling p.value
## Ftest 3.6465 6.0000 127.5621 0.9912 0.002251 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modsininter<-lmer(rel.fitness~damage+pop+(1|pop:fam), REML=F)
p.inter<-KRmodcomp(modtodos,modsininter)
p.inter
## F-test with Kenward-Roger approximation; computing time: 0.07 sec.
## large : rel.fitness ~ damage + pop + (1 | pop:fam) + damage:pop
## small : rel.fitness ~ damage + pop + (1 | pop:fam)
## stat ndf ddf F.scaling p.value
## Ftest 0.6294 3.0000 328.2241 1 0.5965
#Podemos mejorar la presición de la estimación con un boostrap paramétrico
pboost.damage<-PBmodcomp(modtodos,modsindamage) #Demora un minuto o dos dependiendo del equipo
summary (pboost.damage)
## Parametric bootstrap test; time: 20.11 sec; samples: 1000 extremes: 148;
## large : rel.fitness ~ damage + pop + (1 | pop:fam) + damage:pop
## small : rel.fitness ~ pop + (1 | pop:fam)
## stat df ddf p.value
## PBtest 7.0861 0.1489
## Gamma 7.0861 0.1345
## Bartlett 6.8846 4.0000 0.1421
## F 1.7715 4.0000 2.6416 0.3507
## LRT 7.0861 4.0000 0.1314
pboost.pop<-PBmodcomp(modtodos,modsinpop)
summary (pboost.pop)
## Parametric bootstrap test; time: 19.84 sec; samples: 1000 extremes: 3;
## large : rel.fitness ~ damage + pop + (1 | pop:fam) + damage:pop
## small : rel.fitness ~ damage + (1 | pop:fam)
## stat df ddf p.value
## PBtest 20.0241 0.003996 **
## Gamma 20.0241 0.004499 **
## Bartlett 19.0565 6.0000 0.004069 **
## F 3.3373 6.0000 2.377 0.215621
## LRT 20.0241 6.0000 0.002742 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pboost.inter<-PBmodcomp(modtodos,modsininter)
summary (pboost.inter)
## Parametric bootstrap test; time: 21.06 sec; samples: 1000 extremes: 578;
## large : rel.fitness ~ damage + pop + (1 | pop:fam) + damage:pop
## small : rel.fitness ~ damage + pop + (1 | pop:fam)
## stat df ddf p.value
## PBtest 1.9059 0.5784
## Gamma 1.9059 0.5820
## Bartlett 1.9232 3.00000 0.5885
## F 0.6353 3.00000 3.0137 0.6406
## LRT 1.9059 3.00000 0.5922
##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:
## rel.fitness ~ damage + pop + (1 | pop:fam) + damage:pop
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 10 -225.22 470.44
## (1 | pop:fam) 9 -231.55 481.11 12.668 1 0.0003719 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Ajustando valores de p para efectos aleatorios con Restricted Likelihood Ratio Test (diseño no balanceado)
detach("package:lmerTest", unload=TRUE)
library(RLRsim)
#Probando fam
exactRLRT(fit1)
##
## simulated finite sample distribution of RLRT.
##
## (p-value based on 10000 simulated values)
##
## data:
## RLRT = 12.668, p-value < 2.2e-16
#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:fam`
## (Intercept)
## morelos:1 -0.0926386133
## morelos:3 -0.1051733187
## morelos:4 0.0031250953
## morelos:5 0.1375195128
## morelos:7 0.0383546999
## morelos:8 -0.0399333102
## morelos:9 -0.0029374933
## morelos:10 -0.1208309208
## morelos:15 -0.0281215046
## morelos:17 0.1270125072
## morelos:18 0.2341433293
## morelos:19 -0.0111591023
## morelos:25 -0.1694081942
## morelos:26 0.0568560291
## morelos:27 -0.0268087162
## teotihuacan:1 -0.0619715587
## teotihuacan:2 0.1882162457
## teotihuacan:4 -0.0261037200
## teotihuacan:6 -0.1118065941
## teotihuacan:7 0.0824437375
## teotihuacan:8 -0.1520381698
## teotihuacan:9 -0.2294035027
## teotihuacan:11 -0.0743336143
## teotihuacan:12 -0.1148137126
## teotihuacan:13 -0.0910946849
## teotihuacan:16 0.2096156955
## teotihuacan:18 -0.0368229273
## teotihuacan:19 0.1568173472
## teotihuacan:21 -0.0829287516
## teotihuacan:22 -0.0480000968
## teotihuacan:23 0.0617251972
## teotihuacan:24 0.0811389373
## teotihuacan:25 0.0710139095
## teotihuacan:27 0.1985709527
## teotihuacan:29 -0.0202246897
## valdeflores:2 0.0365690987
## valdeflores:3 0.0613603019
## valdeflores:8 -0.0451212866
## valdeflores:9 -0.1630653271
## valdeflores:15 -0.0056896945
## valdeflores:17 0.0913718846
## valdeflores:18 0.1057235077
## valdeflores:24 -0.0286709096
## valdeflores:27 0.0009845376
## valdeflores:28 -0.0534621128
## zubia:5 -0.0508296451
## zubia:12 0.0383132678
## zubia:13 0.0330940670
## zubia:15 -0.0709504939
## zubia:18 -0.1125844785
## zubia:19 0.1595786275
## zubia:22 0.0033786553
dput(cl<-ranef(fit1,condVar=T))
## structure(list(`pop:fam` = structure(list(`(Intercept)` = c(-0.0926386132531279,
## -0.105173318715113, 0.0031250953283876, 0.137519512755068, 0.0383546999254235,
## -0.0399333102241156, -0.00293749333067172, -0.12083092075824,
## -0.0281215046373816, 0.127012507228434, 0.234143329254515, -0.0111591023358657,
## -0.169408194157876, 0.0568560290793499, -0.0268087161587725,
## -0.0619715587086689, 0.188216245689611, -0.0261037199613149,
## -0.111806594133577, 0.0824437374684038, -0.152038169786654, -0.229403502739171,
## -0.0743336143244895, -0.114813712581708, -0.0910946848678053,
## 0.209615695474272, -0.0368229273312455, 0.156817347173257, -0.0829287515821665,
## -0.0480000968401863, 0.0617251972399262, 0.0811389373422072,
## 0.0710139095457095, 0.19857095266566, -0.0202246897420698, 0.0365690987148599,
## 0.0613603019057468, -0.0451212865643001, -0.163065327070625,
## -0.00568969454657004, 0.0913718846410781, 0.105723507716639,
## -0.0286709096439452, 0.000984537615846747, -0.053462112768736,
## -0.050829645077446, 0.038313267782879, 0.0330940670137085, -0.0709504939206507,
## -0.112584478539827, 0.159578627475502, 0.00337865526582859)), class = "data.frame", row.names = c("morelos:1",
## "morelos:3", "morelos:4", "morelos:5", "morelos:7", "morelos:8",
## "morelos:9", "morelos:10", "morelos:15", "morelos:17", "morelos:18",
## "morelos:19", "morelos:25", "morelos:26", "morelos:27", "teotihuacan:1",
## "teotihuacan:2", "teotihuacan:4", "teotihuacan:6", "teotihuacan:7",
## "teotihuacan:8", "teotihuacan:9", "teotihuacan:11", "teotihuacan:12",
## "teotihuacan:13", "teotihuacan:16", "teotihuacan:18", "teotihuacan:19",
## "teotihuacan:21", "teotihuacan:22", "teotihuacan:23", "teotihuacan:24",
## "teotihuacan:25", "teotihuacan:27", "teotihuacan:29", "valdeflores:2",
## "valdeflores:3", "valdeflores:8", "valdeflores:9", "valdeflores:15",
## "valdeflores:17", "valdeflores:18", "valdeflores:24", "valdeflores:27",
## "valdeflores:28", "zubia:5", "zubia:12", "zubia:13", "zubia:15",
## "zubia:18", "zubia:19", "zubia:22"), postVar = structure(c(0.0118566658109377,
## 0.0110729855345261, 0.0110729855345261, 0.0110729855345261, 0.0110729855345261,
## 0.0110729855345261, 0.0127597233942188, 0.0110729855345261, 0.0110729855345261,
## 0.0110729855345261, 0.0110729855345261, 0.0110729855345261, 0.0118566658109377,
## 0.0118566658109377, 0.0103864789860161, 0.0118566658109377, 0.0118566658109377,
## 0.0127597233942188, 0.0127597233942188, 0.0110729855345261, 0.0118566658109377,
## 0.0110729855345261, 0.0150526853182835, 0.0127597233942188, 0.0110729855345261,
## 0.0110729855345261, 0.0118566658109377, 0.0110729855345261, 0.0118566658109377,
## 0.0118566658109377, 0.0118566658109377, 0.0110729855345261, 0.0138116840570645,
## 0.0110729855345261, 0.0165387144996812, 0.0110729855345261, 0.0110729855345261,
## 0.0118566658109377, 0.0110729855345261, 0.0110729855345261, 0.0150526853182835,
## 0.0110729855345261, 0.0110729855345261, 0.0118566658109377, 0.0110729855345261,
## 0.0118566658109377, 0.0110729855345261, 0.0110729855345261, 0.0118566658109377,
## 0.0118566658109377, 0.0118566658109377, 0.0118566658109377), .Dim = c(1L,
## 1L, 52L)))), class = "ranef.mer")
dotplot(cl)
## $`pop:fam`
##Podemos hacer analisis a posteriori sobre pop ya que es fijo (sin basarse en estimación como cuando es aleatorio)
library(emmeans)
## NOTE: As of emmeans versions > 1.2.3,
## The 'cld' function will be deprecated in favor of 'CLD'.
## You may use 'cld' only if you have package:multcomp attached.
multc<-emmeans(fit1, "pop")
## NOTE: Results may be misleading due to involvement in interactions
CLD(multc)
## pop emmean SE df lower.CL upper.CL .group
## valdeflores 0.8068450 0.06827641 46.41 0.6694446 0.9442454 1
## zubia 0.8084817 0.08161305 47.32 0.6443268 0.9726366 1
## teotihuacan 1.0902301 0.04949436 50.58 0.9908461 1.1896141 2
## morelos 1.1020497 0.05494535 44.54 0.9913526 1.2127468 2
##
## Results are averaged over the levels of: damage
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
plot(multc, comparisons = T)
