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.fit.origen~damage+(1|pop)+(1|pop:fam))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: rel.fit.origen ~ damage + (1 | pop) + (1 | pop:fam)
##
## REML criterion at convergence: 407.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4975 -0.4383 0.0682 0.4839 4.8566
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 0.02069 0.1439
## pop (Intercept) 0.00000 0.0000
## Residual 0.15440 0.3929
## Number of obs: 376, groups: pop:fam, 52; pop, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.05174 0.03512 29.946
## damage1 -0.10296 0.04079 -2.524
##
## Correlation of Fixed Effects:
## (Intr)
## damage1 -0.583
##Probando efectos fijos (damage) con la correción Kenward-Roger para datos desbalanceados
library(pbkrtest)
modtodos<-lmer(rel.fit.origen~damage+(1|pop)+(1|pop:fam), REML = F)
modsindamage<-lmer(rel.fit.origen~1+(1|pop)+(1|pop:fam), REML = F)
p.damage<-KRmodcomp(modtodos,modsindamage)
p.damage
## F-test with Kenward-Roger approximation; computing time: 0.17 sec.
## large : rel.fit.origen ~ damage + (1 | pop) + (1 | pop:fam)
## small : rel.fit.origen ~ 1 + (1 | pop) + (1 | pop:fam)
## stat ndf ddf F.scaling p.value
## Ftest 6.3568 1.0000 332.4978 1 0.01216 *
## ---
## 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: 23.81 sec; samples: 1000 extremes: 14;
## large : rel.fit.origen ~ damage + (1 | pop) + (1 | pop:fam)
## small : rel.fit.origen ~ 1 + (1 | pop) + (1 | pop:fam)
## stat df ddf p.value
## PBtest 6.3111 0.01499 *
## Gamma 6.3111 0.01449 *
## Bartlett 6.0349 1.0000 0.01403 *
## F 6.3111 1.0000 45.701 0.01559 *
## LRT 6.3111 1.0000 0.01200 *
## ---
## 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.fit.origen ~ damage + (1 | pop) + (1 | pop:fam)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 5 -203.75 417.50
## (1 | pop) 4 -203.75 415.50 0.000 1 0.9999997
## (1 | pop:fam) 4 -209.91 427.81 12.314 1 0.0004496 ***
## ---
## 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.fit.origen~damage+(1|pop:fam))
mod.solo.pop<-lmer(rel.fit.origen~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 = 1.1369e-13, p-value = 0.3752
#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.314, p-value = 3e-04
#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.080502061
## morelos:3 -0.092532154
## morelos:4 0.003853304
## morelos:5 0.123464174
## morelos:7 0.035207610
## morelos:8 -0.034468625
## morelos:9 -0.001682416
## morelos:10 -0.106467401
## morelos:15 -0.025486917
## morelos:17 0.115643737
## morelos:18 0.207928470
## morelos:19 -0.010390399
## morelos:25 -0.150359917
## morelos:26 0.052346215
## morelos:27 -0.023477700
## teotihuacan:1 -0.061878471
## teotihuacan:2 0.170877619
## teotihuacan:4 -0.027024722
## teotihuacan:6 -0.103051641
## teotihuacan:7 0.068884457
## teotihuacan:8 -0.136702618
## teotihuacan:9 -0.208659212
## teotihuacan:11 -0.068679222
## teotihuacan:12 -0.105719253
## teotihuacan:13 -0.085564530
## teotihuacan:16 0.177191077
## teotihuacan:18 -0.039530163
## teotihuacan:19 0.135076880
## teotihuacan:21 -0.080502061
## teotihuacan:22 -0.049462744
## teotihuacan:23 0.053258001
## teotihuacan:24 0.067723186
## teotihuacan:25 0.062411572
## teotihuacan:27 0.172237538
## teotihuacan:29 -0.009489636
## valdeflores:2 0.046650090
## valdeflores:3 0.076847099
## valdeflores:8 -0.052738139
## valdeflores:9 -0.196515299
## valdeflores:15 -0.004509932
## valdeflores:17 0.111865133
## valdeflores:18 0.130883853
## valdeflores:24 -0.032815723
## valdeflores:27 0.003336049
## valdeflores:28 -0.063012733
## zubia:5 -0.044564828
## zubia:12 0.049828723
## zubia:13 0.043471458
## zubia:15 -0.088098934
## zubia:18 -0.129202940
## zubia:19 0.192272004
## zubia:22 0.011832138
##
## $pop
## (Intercept)
## morelos 0
## teotihuacan 0
## valdeflores 0
## zubia 0
dput(cl<-ranef(fit1,condVar=T))
## structure(list(`pop:fam` = structure(list(`(Intercept)` = c(-0.0805020607974984,
## -0.0925321537094571, 0.00385330437305278, 0.123464174083789,
## 0.035207610116433, -0.0344686247682464, -0.00168241607126479,
## -0.106467400634651, -0.0254869167238369, 0.115643736627259, 0.207928469822918,
## -0.0103903992107641, -0.150359916943788, 0.0523462150249431,
## -0.0234776999950628, -0.0618784707909731, 0.170877618775708,
## -0.0270247223408065, -0.103051641372293, 0.0688844567768601,
## -0.136702617522689, -0.208659211591878, -0.0686792216798531,
## -0.105719252582019, -0.0855645302468602, 0.177191077267429, -0.0395301626033523,
## 0.135076879897902, -0.0805020609357988, -0.049462743981656, 0.0532580012768604,
## 0.067723186253659, 0.0624115721656215, 0.172237538451322, -0.00948963581795184,
## 0.0466500902452007, 0.0768470993293142, -0.0527381385483701,
## -0.196515299301221, -0.00450993233028608, 0.111865133312319,
## 0.13088385252592, -0.0328157232906337, 0.00333604909234633, -0.0630127325041027,
## -0.0445648278009957, 0.0498287226339945, 0.0434714575976958,
## -0.0880989340972051, -0.129202939807469, 0.192272004175527, 0.0118321381749186
## )), 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.0106771175731766, 0.00998654592798574,
## 0.00998654592798574, 0.00998654592798574, 0.00998654592798574,
## 0.00998654592798574, 0.0114702904482406, 0.00998654592798574,
## 0.00998654592798574, 0.00998654592798574, 0.00998654592798574,
## 0.00998654592798574, 0.0106771175731766, 0.0106771175731766,
## 0.00937987685753643, 0.0106771175731766, 0.0106771175731766,
## 0.0114702904482406, 0.0114702904482406, 0.00998654592798574,
## 0.0106771175731766, 0.00998654592798574, 0.0134718637865095,
## 0.0114702904482406, 0.00998654592798574, 0.00998654592798574,
## 0.0106771175731766, 0.00998654592798574, 0.0106771175731766,
## 0.0106771175731766, 0.0106771175731766, 0.00998654592798574,
## 0.0123907653730333, 0.00998654592798574, 0.0147596479996296,
## 0.00998654592798574, 0.00998654592798574, 0.0106771175731766,
## 0.00998654592798574, 0.00998654592798574, 0.0134718637865095,
## 0.00998654592798574, 0.00998654592798574, 0.0106771175731766,
## 0.00998654592798574, 0.0106771175731766, 0.00998654592798574,
## 0.00998654592798574, 0.0106771175731766, 0.0106771175731766,
## 0.0106771175731766, 0.0106771175731766), .Dim = c(1L, 1L, 52L
## ))), pop = structure(list(`(Intercept)` = c(0, 0, 0, 0)), class = "data.frame", row.names = c("morelos",
## "teotihuacan", "valdeflores", "zubia"), postVar = structure(c(0,
## 0, 0, 0), .Dim = c(1L, 1L, 4L)))), class = "ranef.mer")
dotplot(cl)
## $`pop:fam`
##
## $pop
