Cargamos paquetes:

library(lme4)
library(nlme)
library(MCMCglmm)
library(AICcmodavg)
library(mgcv)
library(MuMIn)
library(MASS)
library(lattice)
library(lavaan)
library(piecewiseSEM)
library(car)

Cargamos la base de datos:

load("res2.Rda")

Estructura de variables:

str(res2)
## Classes 'data.table' and 'data.frame':   272 obs. of  7 variables:
##  $ Site     : Factor w/ 36 levels "CE01","CE02",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ mean     : num  1.31e+03 2.76e+03 6.19e+02 3.83e+01 6.11e-02 ...
##  $ sd       : num  1101.67 1310.07 1226.58 98.61 0.32 ...
##  $ Frec_Dom : chr  "F12" "F23" "F34" "F45" ...
##  $ prom.vert: num  0.142 0.142 0.142 0.142 0.142 ...
##  $ prom.hor : num  0.0542 0.0542 0.0542 0.0542 0.0542 ...
##  $ AB       : num  20296 20296 20296 20296 20296 ...
##  - attr(*, ".internal.selfref")=<externalptr> 
##  - attr(*, "sorted")= chr "Site"

Estandarización de variables:

Site <- as.factor(res2$Site)
FDm <- as.numeric(res2$mean)
FD <- as.factor(res2$Frec_Dom)
Fv <- as.numeric(((res2$prom.vert)-mean(res2$prom.vert))/sd(res2$prom.vert))
Fh <- as.numeric(((res2$prom.hor)-mean(res2$prom.hor))/sd(res2$prom.hor))
AB <- as.numeric(((res2$AB)-mean(res2$AB))/sd(res2$AB))

Agrupación de variables transformadas:

data<- data.frame(Site,FDm, FD, Fv, Fh,AB )
attach(data)
## The following objects are masked _by_ .GlobalEnv:
## 
##     AB, FD, FDm, Fh, Fv, Site
mod1 <- lm(FDm~Fv+Fh+AB,data = res2)
mod2 <- lm(FDm~Fh,data = res2)
mod3 <- lm(FDm~Fv,data = res2)
mod4 <- lm(FDm~Fv+Fh+Fv:Fh ,data = res2)
mod5 <- lm(FDm~AB)
mod0 <- lm(FDm~1,data = res2)
AIC(mod1, mod2, mod3, mod4, mod5)
##      df      AIC
## mod1  5 4465.891
## mod2  3 4463.353
## mod3  3 4462.030
## mod4  5 4465.943
## mod5  3 4463.010
AIC(mod1, mod0)
##      df      AIC
## mod1  5 4465.891
## mod0  2 4461.504

No hay fuertes evidencias de que existe un mejor modelo.

#Comprobación de Supuestos
par(mfrow=c(2,2))
plot(mod2)

#Revizamos la ortogonalidad de las variables independientes

cor_PCA <- (princomp(res2[,5:6]))
biplot(cor_PCA)

#Seleccion de variables a traves del método “backward”

step <- stepAIC(mod1, direction="backward")
## Start:  AIC=3691.99
## FDm ~ Fv + Fh + AB
## 
##        Df Sum of Sq       RSS    AIC
## - Fh    1     29473 207371016 3690.0
## - AB    1     64042 207405584 3690.1
## - Fv    1    794472 208136015 3691.0
## <none>              207341542 3692.0
## 
## Step:  AIC=3690.03
## FDm ~ Fv + AB
## 
##        Df Sum of Sq       RSS    AIC
## - AB    1     76441 207447457 3688.1
## - Fv    1    825430 208196445 3689.1
## <none>              207371016 3690.0
## 
## Step:  AIC=3688.13
## FDm ~ Fv
## 
##        Df Sum of Sq       RSS    AIC
## - Fv    1   1127158 208574615 3687.6
## <none>              207447457 3688.1
## 
## Step:  AIC=3687.6
## FDm ~ 1

No se muetsra una efecto importante para eliminar alguna variable

#Para la importancia relativa de cada variable del modelo

library(leaps)
attach(data)
## The following objects are masked _by_ .GlobalEnv:
## 
##     AB, FD, FDm, Fh, Fv, Site
## The following objects are masked from data (pos = 4):
## 
##     AB, FD, FDm, Fh, Fv, Site
leaps <- regsubsets(FDm~Fv+Fh,data = data, nbest = 10)
summary(leaps)
## Subset selection object
## Call: regsubsets.formula(FDm ~ Fv + Fh, data = data, nbest = 10)
## 2 Variables  (and intercept)
##    Forced in Forced out
## Fv     FALSE      FALSE
## Fh     FALSE      FALSE
## 10 subsets of each size up to 2
## Selection Algorithm: exhaustive
##          Fv  Fh 
## 1  ( 1 ) "*" " "
## 1  ( 2 ) " " "*"
## 2  ( 1 ) "*" "*"
plot(leaps, scale = "r2")

#intercepto y Fh explica 0.0056 de la varianza

#Autocorrelación

vif(mod1) #Si es menor a 2-4 esta bien
##       Fv       Fh       AB 
## 1.133011 1.032976 1.142911
sqrt(vif(mod1)) #Muetras la inflacion de de cada variable (entre 2 a 4)
##       Fv       Fh       AB 
## 1.064430 1.016354 1.069070
max(vif(mod1)) #Si es menor a 2-4 esta bien
## [1] 1.142911

#Relación entre variables

par(mfrow=c(1,3))
library(visreg)
visreg(mod1)

#Modelo lineal genralizado
glm1 <- glm((FDm+1)~Fv+Fh+AB, family=poisson, data=data)
summary(glm1)
## 
## Call:
## glm(formula = (FDm + 1) ~ Fv + Fh + AB, family = poisson, data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -34.838  -30.314  -27.133    2.557  100.298  
## 
## Coefficients:
##              Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)  6.108153   0.002877 2122.928  < 2e-16 ***
## Fv          -0.136771   0.003238  -42.242  < 2e-16 ***
## Fh          -0.022443   0.003102   -7.235 4.67e-13 ***
## AB           0.031009   0.002813   11.024  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 319740  on 271  degrees of freedom
## Residual deviance: 316954  on 268  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 6
glm1quasi<- glm((FDm+1)~Fv+Fh+AB, family=quasipoisson, data) #Modelo de quasi, por la sobredispersión
summary(glm1quasi)
## 
## Call:
## glm(formula = (FDm + 1) ~ Fv + Fh + AB, family = quasipoisson, 
##     data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -34.838  -30.314  -27.133    2.557  100.298  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.10815    0.11767  51.910   <2e-16 ***
## Fv          -0.13677    0.13241  -1.033    0.303    
## Fh          -0.02244    0.12687  -0.177    0.860    
## AB           0.03101    0.11503   0.270    0.788    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1672.515)
## 
##     Null deviance: 319740  on 271  degrees of freedom
## Residual deviance: 316954  on 268  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 6

#Modelos Mixtos

M1 <- lmer(log(FDm+1)~Fv+poly(AB,2)+Fh+poly(AB,2)+(1|Site),data=data)
M2 <- lmer(log(FDm+1)~Fv+poly(AB,2)+(1|Site),data=data)
M3 <- lmer(log(FDm+1)~Fh+poly(AB,2)+(1|Site),data=data)
M0 <- lmer(log(FDm+1)~1+(1|Site),data=data)
AIC(M1,M2,M3)
##    df      AIC
## M1  7 1390.601
## M2  6 1387.175
## M3  6 1387.222

#Elección del Modelo

M1
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(FDm + 1) ~ Fv + poly(AB, 2) + Fh + poly(AB, 2) + (1 | Site)
##    Data: data
## REML criterion at convergence: 1376.601
## Random effects:
##  Groups   Name        Std.Dev.
##  Site     (Intercept) 0.000   
##  Residual             3.089   
## Number of obs: 272, groups:  Site, 34
## Fixed Effects:
##  (Intercept)            Fv  poly(AB, 2)1  poly(AB, 2)2            Fh  
##     3.013469     -0.007105      1.419687     -0.985677      0.043415
### Goodness of fitness ####
library(MuMIn)
r.squaredGLMM(M1)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help
## page.
##              R2m         R2c
## [1,] 0.001257965 0.001257965
plot(M1)