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.pop~damage+(1|pop)+(1|pop:fam))
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: rel.fit.pop ~ damage + (1 | pop) + (1 | pop:fam)
##
## REML criterion at convergence: 407.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4950 -0.4398 0.0655 0.4854 4.8635
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:fam (Intercept) 0.02071 0.1439
## pop (Intercept) 0.00000 0.0000
## Residual 0.15447 0.3930
## Number of obs: 376, groups: pop:fam, 52; pop, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.05175 0.03513 29.937
## damage1 -0.10296 0.04080 -2.523
##
## 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.pop~damage+(1|pop)+(1|pop:fam), REML = F)
modsindamage<-lmer(rel.fit.pop~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.fit.pop ~ damage + (1 | pop) + (1 | pop:fam)
## small : rel.fit.pop ~ 1 + (1 | pop) + (1 | pop:fam)
## stat ndf ddf F.scaling p.value
## Ftest 6.3537 1.0000 332.4956 1 0.01218 *
## ---
## 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
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge: degenerate Hessian with 1
## negative eigenvalues
summary (pboost.damage)
## Parametric bootstrap test; time: 23.89 sec; samples: 1000 extremes: 18;
## large : rel.fit.pop ~ damage + (1 | pop) + (1 | pop:fam)
## small : rel.fit.pop ~ 1 + (1 | pop) + (1 | pop:fam)
## stat df ddf p.value
## PBtest 6.3081 0.01898 *
## Gamma 6.3081 0.01454 *
## Bartlett 5.9483 1.0000 0.01473 *
## F 6.3081 1.0000 35.066 0.01678 *
## LRT 6.3081 1.0000 0.01202 *
## ---
## 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.pop ~ damage + (1 | pop) + (1 | pop:fam)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 5 -203.83 417.67
## (1 | pop) 4 -203.83 415.67 0.000 1 0.9999998
## (1 | pop:fam) 4 -210.00 428.00 12.331 1 0.0004454 ***
## ---
## 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.pop~damage+(1|pop:fam))
mod.solo.pop<-lmer(rel.fit.pop~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 = 5.6843e-14, p-value = 0.381
#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.331, p-value = 2e-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.081177821
## morelos:3 -0.093238404
## morelos:4 0.003020803
## morelos:5 0.122474999
## morelos:7 0.034334039
## morelos:8 -0.035250930
## morelos:9 -0.002393108
## morelos:10 -0.107155398
## morelos:15 -0.026270644
## morelos:17 0.114654463
## morelos:18 0.206839001
## morelos:19 -0.011193901
## morelos:25 -0.150934480
## morelos:26 0.051499039
## morelos:27 -0.024316008
## teotihuacan:1 -0.061323868
## teotihuacan:2 0.171797662
## teotihuacan:4 -0.026469612
## teotihuacan:6 -0.102624086
## teotihuacan:7 0.069694137
## teotihuacan:8 -0.136281781
## teotihuacan:9 -0.208303494
## teotihuacan:11 -0.068328171
## teotihuacan:12 -0.105296173
## teotihuacan:13 -0.085007473
## teotihuacan:16 0.178187180
## teotihuacan:18 -0.038938568
## teotihuacan:19 0.135994827
## teotihuacan:21 -0.079978285
## teotihuacan:22 -0.048887590
## teotihuacan:23 0.053993269
## teotihuacan:24 0.068530967
## teotihuacan:25 0.063052706
## teotihuacan:27 0.173216267
## teotihuacan:29 -0.009140644
## valdeflores:2 0.045622611
## valdeflores:3 0.075772631
## valdeflores:8 -0.053556831
## valdeflores:9 -0.197164394
## valdeflores:15 -0.005469784
## valdeflores:17 0.111057958
## valdeflores:18 0.129725299
## valdeflores:24 -0.033719547
## valdeflores:27 0.002431197
## valdeflores:28 -0.063869567
## zubia:5 -0.043440943
## zubia:12 0.051275936
## zubia:13 0.044900780
## zubia:15 -0.087135027
## zubia:18 -0.128337210
## zubia:19 0.194030435
## zubia:22 0.013097537
##
## $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.081177820602813,
## -0.0932384039557863, 0.00302080277412847, 0.122474999349552,
## 0.0343340388682421, -0.0352509300213233, -0.00239310751772336,
## -0.107155397811335, -0.0262706443377932, 0.114654463213289, 0.20683900120611,
## -0.0111939011937619, -0.150934480473795, 0.0514990390632221,
## -0.024316008488634, -0.061323868013736, 0.171797661997535, -0.0264696123868726,
## -0.102624085955583, 0.0696941368313197, -0.136281780504614, -0.20830349352032,
## -0.068328170786223, -0.105296172802223, -0.0850074733410827,
## 0.178187180303525, -0.038938568004839, 0.13599482696066, -0.0799782846186458,
## -0.0488875901387873, 0.0539932690868138, 0.0685309667405191,
## 0.0630527058196984, 0.173216267019652, -0.00914064421483926,
## 0.0456226111454052, 0.0757726314086286, -0.0535568307853157,
## -0.197164393839297, -0.00546978426957716, 0.111057958230407,
## 0.129725299125499, -0.0337195472922027, 0.0024311966093401, -0.0638695674907302,
## -0.0434409433164338, 0.0512759356744959, 0.0449007799333705,
## -0.0871350274277879, -0.128337210301728, 0.194030434737237, 0.013097537315193
## )), 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.0106847693537415, 0.00999350724937945,
## 0.00999350724937945, 0.00999350724937945, 0.00999350724937945,
## 0.00999350724937945, 0.0114787686553023, 0.00999350724937945,
## 0.00999350724937945, 0.00999350724937945, 0.00999350724937945,
## 0.00999350724937945, 0.0106847693537415, 0.0106847693537415,
## 0.00938625393115459, 0.0106847693537415, 0.0106847693537415,
## 0.0114787686553023, 0.0114787686553023, 0.00999350724937945,
## 0.0106847693537415, 0.00999350724937945, 0.0134825861736159,
## 0.0114787686553023, 0.00999350724937945, 0.00999350724937945,
## 0.0106847693537415, 0.00999350724937945, 0.0106847693537415,
## 0.0106847693537415, 0.0106847693537415, 0.00999350724937945,
## 0.0124002473922464, 0.00999350724937945, 0.0147719344437405,
## 0.00999350724937945, 0.00999350724937945, 0.0106847693537415,
## 0.00999350724937945, 0.00999350724937945, 0.0134825861736159,
## 0.00999350724937945, 0.00999350724937945, 0.0106847693537415,
## 0.00999350724937945, 0.0106847693537415, 0.00999350724937945,
## 0.00999350724937945, 0.0106847693537415, 0.0106847693537415,
## 0.0106847693537415, 0.0106847693537415), .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
