por.origen.R

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*origen+(1|pop))
summary(fit1) 

## Linear mixed model fit by REML ['lmerMod']
## Formula: rel.fit.origen ~ damage * origen + (1 | pop)
## 
## REML criterion at convergence: 426.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4024 -0.5012  0.0562  0.5698  4.9417 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  pop      (Intercept) 0.0000   0.0000  
##  Residual             0.1755   0.4189  
## Number of obs: 376, groups:  pop, 4
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)        1.08009    0.05156  20.947
## damage1           -0.15423    0.07505  -2.055
## origenMEX         -0.04431    0.06401  -0.692
## damage1:origenMEX  0.08272    0.09182   0.901
## 
## Correlation of Fixed Effects:
##             (Intr) damag1 orgMEX
## damage1     -0.687              
## origenMEX   -0.806  0.553       
## dmg1:rgnMEX  0.562 -0.817 -0.697

##Probando efectos fijos (damage) con la correción Kenward-Roger para datos desbalanceados
library(pbkrtest)
modtodos<-lmer(rel.fit.origen~damage+origen+(1|pop)+damage:origen, REML = F)
modsindamage<-lmer(rel.fit.origen~origen+(1|pop), REML = F)
p.damage<-KRmodcomp(modtodos,modsindamage)
p.damage

## F-test with Kenward-Roger approximation; computing time: 0.13 sec.
## large : rel.fit.origen ~ damage + origen + (1 | pop) + damage:origen
## small : rel.fit.origen ~ origen + (1 | pop)
##           stat      ndf      ddf F.scaling p.value  
## Ftest   3.0251   2.0000 370.2336         1 0.04976 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

modsinorigen<-lmer(rel.fit.origen~damage+(1|pop), REML = F)
p.origen<-KRmodcomp(modtodos,modsinorigen)
p.origen

## F-test with Kenward-Roger approximation; computing time: 0.08 sec.
## large : rel.fit.origen ~ damage + origen + (1 | pop) + damage:origen
## small : rel.fit.origen ~ damage + (1 | pop)
##         stat    ndf    ddf F.scaling p.value
## Ftest 0.3715 2.0000 6.3560   0.90687  0.7038

modsininter<-lmer(rel.fit.origen~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 : rel.fit.origen ~ damage + origen + (1 | pop) + damage:origen
## small : rel.fit.origen ~ damage + origen + (1 | pop)
##           stat      ndf      ddf F.scaling p.value
## Ftest   0.8115   1.0000 370.2284         1  0.3683

##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 + origen + (1 | pop) + damage:origen
##           npar  logLik    AIC LRT Df Pr(>Chisq)
## <none>       6 -213.13 438.26                  
## (1 | pop)    5 -213.13 436.26   0  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               0
## teotihuacan           0
## valdeflores           0
## zubia                 0

dput(cl<-ranef(fit1,condVar=T))

## structure(list(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