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+(1|pop)+(1|pop:fam))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: rel.fitness ~ damage + (1 | pop) + (1 | pop:fam)
##
## REML criterion at convergence: 443.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5304 -0.4240 0.0343 0.4841 5.1236
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 0.02362 0.1537
## pop (Intercept) 0.02335 0.1528
## Residual 0.16693 0.4086
## Number of obs: 376, groups: pop:fam, 52; pop, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.00829 0.08546 11.799
## damage1 -0.09631 0.04244 -2.269
##
## Correlation of Fixed Effects:
## (Intr)
## damage1 -0.246
##Probando efectos fijos (damage) con la correción Kenward-Roger para datos desbalanceados
library(pbkrtest)
modtodos<-lmer(rel.fitness~damage+(1|pop)+(1|pop:fam), REML = F)
modsindamage<-lmer(rel.fitness~1+(1|pop)+(1|pop:fam), REML = F)
p.damage<-KRmodcomp(modtodos,modsindamage)
p.damage
## F-test with Kenward-Roger approximation; computing time: 0.16 sec.
## large : rel.fitness ~ damage + (1 | pop) + (1 | pop:fam)
## small : rel.fitness ~ 1 + (1 | pop) + (1 | pop:fam)
## stat ndf ddf F.scaling p.value
## Ftest 5.1432 1.0000 331.8785 1 0.02398 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#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: 24.89 sec; samples: 1000 extremes: 19;
## large : rel.fitness ~ damage + (1 | pop) + (1 | pop:fam)
## small : rel.fitness ~ 1 + (1 | pop) + (1 | pop:fam)
## stat df ddf p.value
## PBtest 5.0879 0.01998 *
## Gamma 5.0879 0.02249 *
## Bartlett 5.1392 1.0000 0.02339 *
## F 5.0879 1.0000 -198.47
## LRT 5.0879 1.0000 0.02409 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##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 + (1 | pop) + (1 | pop:fam)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 5 -221.92 453.84
## (1 | pop) 4 -226.35 460.70 8.860 1 0.0029149 **
## (1 | pop:fam) 4 -228.35 464.70 12.854 1 0.0003368 ***
## ---
## 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)
mod.solo.fam<-lmer(rel.fitness~damage+(1|pop:fam))
mod.solo.pop<-lmer(rel.fitness~damage+(1|pop))
library(RLRsim)
#Probando pop
exactRLRT(mod.solo.pop, fit1, mod.solo.fam)
##
## simulated finite sample distribution of RLRT.
##
## (p-value based on 10000 simulated values)
##
## data:
## RLRT = 8.86, p-value = 4e-04
#Probando fam
exactRLRT(mod.solo.fam, fit1, mod.solo.pop)
##
## simulated finite sample distribution of RLRT.
##
## (p-value based on 10000 simulated values)
##
## data:
## RLRT = 12.854, 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.083121102
## morelos:3 -0.096642723
## morelos:4 0.012102787
## morelos:5 0.147052035
## morelos:7 0.047477833
## morelos:8 -0.031133380
## morelos:9 0.004801191
## morelos:10 -0.112364966
## morelos:15 -0.022121008
## morelos:17 0.139349851
## morelos:18 0.241226553
## morelos:19 -0.005088579
## morelos:25 -0.163280814
## morelos:26 0.067034583
## morelos:27 -0.018790278
## teotihuacan:1 -0.058927529
## teotihuacan:2 0.201750294
## teotihuacan:4 -0.021001604
## teotihuacan:6 -0.107112444
## teotihuacan:7 0.088827612
## teotihuacan:8 -0.144700538
## teotihuacan:9 -0.224307051
## teotihuacan:11 -0.070637504
## teotihuacan:12 -0.110133877
## teotihuacan:13 -0.085427242
## teotihuacan:16 0.212144115
## teotihuacan:18 -0.033667694
## teotihuacan:19 0.163508264
## teotihuacan:21 -0.079977392
## teotihuacan:22 -0.044894288
## teotihuacan:23 0.070008058
## teotihuacan:24 0.087517425
## teotihuacan:25 0.078830940
## teotihuacan:27 0.205434244
## teotihuacan:29 -0.007123643
## valdeflores:2 0.023293487
## valdeflores:3 0.048187038
## valdeflores:8 -0.058249092
## valdeflores:9 -0.177165104
## valdeflores:15 -0.019783743
## valdeflores:17 0.082745787
## valdeflores:18 0.092733391
## valdeflores:24 -0.042215856
## valdeflores:27 -0.011939395
## valdeflores:28 -0.067109407
## zubia:5 -0.058561156
## zubia:12 0.019409819
## zubia:13 0.014169072
## zubia:15 -0.092581406
## zubia:18 -0.127494275
## zubia:19 0.138967080
## zubia:22 -0.011018370
##
## $pop
## (Intercept)
## morelos 0.1250333
## teotihuacan 0.1187156
## valdeflores -0.1279993
## zubia -0.1157496
dput(cl<-ranef(fit1,condVar=T))
## structure(list(`pop:fam` = structure(list(`(Intercept)` = c(-0.0831211015961688,
## -0.0966427232367505, 0.0121027872772155, 0.147052035263627, 0.0474778328948827,
## -0.0311333795224529, 0.00480119051608715, -0.112364965733492,
## -0.0221210080504849, 0.139349850555772, 0.241226553438081, -0.00508857859052335,
## -0.163280813567357, 0.0670345834631499, -0.018790278323802, -0.0589275292504622,
## 0.2017502944623, -0.0210016038544219, -0.10711244429896, 0.0888276119938979,
## -0.144700537960365, -0.224307050579486, -0.0706375042252001,
## -0.110133877353403, -0.0854272420352587, 0.212144114714638, -0.0336676944238539,
## 0.163508263853417, -0.0799773916652083, -0.0448942876880227,
## 0.0700080583856966, 0.0875174251634153, 0.0788309399868052, 0.205434243888972,
## -0.0071236434101355, 0.0232934871343225, 0.0481870377099163,
## -0.0582490921736282, -0.177165104433648, -0.0197837429920922,
## 0.0827457873898265, 0.0927333914728053, -0.0422158564472369,
## -0.0119393950033609, -0.0671094069564617, -0.0585611563630049,
## 0.0194098191001365, 0.0141690715127308, -0.0925814059590784,
## -0.127494275070024, 0.138967080389867, -0.0110183698032119)), 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.0125289848439243,
## 0.0118331001164509, 0.0118331001164509, 0.0118331001164509, 0.0118331001164509,
## 0.0118331001164509, 0.0133387581514417, 0.0118331001164509, 0.0118331001164509,
## 0.0118331001164509, 0.0118331001164509, 0.0118331001164509, 0.0125289848439243,
## 0.0125289848439243, 0.0112287185125231, 0.012415378014358, 0.012415378014358,
## 0.0132420282373634, 0.0132420282373634, 0.011703758341071, 0.012415378014358,
## 0.011703758341071, 0.0153729729916582, 0.0132420282373634, 0.011703758341071,
## 0.011703758341071, 0.012415378014358, 0.011703758341071, 0.012415378014358,
## 0.012415378014358, 0.012415378014358, 0.011703758341071, 0.0142139523290543,
## 0.011703758341071, 0.0167785974789612, 0.0121735740850801, 0.0121735740850801,
## 0.0128280387885001, 0.0121735740850801, 0.0121735740850801, 0.0155906716802954,
## 0.0121735740850801, 0.0121735740850801, 0.0128280387885001, 0.0121735740850801,
## 0.0131490219704668, 0.0125390145676741, 0.0125390145676741, 0.0131490219704668,
## 0.0131490219704668, 0.0131490219704668, 0.0131490219704668), .Dim = c(1L,
## 1L, 52L))), pop = structure(list(`(Intercept)` = c(0.125033272547878,
## 0.118715643938804, -0.127999340926246, -0.115749575560472)), class = "data.frame", row.names = c("morelos",
## "teotihuacan", "valdeflores", "zubia"), postVar = structure(c(0.00267402705155347,
## 0.00221522252976355, 0.00388176520034604, 0.00517806521916659
## ), .Dim = c(1L, 1L, 4L)))), class = "ranef.mer")
dotplot(cl)
## $`pop:fam`
##
## $pop
