R linear mixed effects- finding individual fixed effects variances (in %)-CodePudding

(struggling) R user here. I've made some lme4 models, and I am trying to find individual fixed effects variances. I am estimating the variance of the fixed effect components by multiplying the design matrix of the fixed effects with the vector of fixed effect estimates, followed by calculating the variance of these fitted values (As per Nakagawa & Schielzeth, 2013).

I am getting values, however they are sometimes negtive. Is this incorrect?

I also would like to use these values to get the percentage of variance explained, how could this be done? I'm lost when trying to understand all of the maths, and I feel *100 is probably wrong...

Here is an example of my code:

data
     Ho          SP        SL          DP         DL            species
1 0.432 -0.62240975 0.9615917 -0.20535833 -0.8038894 Accipiter gentilis
2 0.355 -0.85728592 1.0821313 -0.07917467 -0.8055025 Accipiter gentilis
3 0.370 -0.91221808 0.7616323 -0.17403906 -0.8041293 Accipiter gentilis
4 0.447 -1.15202246 0.9326085 -0.07640844 -0.8031998 Accipiter gentilis
5 0.452  0.12952583 0.7566257  0.35393076 -0.8058789 Ales Alces
6 0.406 -0.06462976 0.6679433  0.86746704 -0.7952550 Ales Alces

lme <- lmer(data$Ho ~ data$SP   data$DP   data$SL   data$DL   (1|data$spp)) 
vector <- fixef(lme)
vector <- vector[-1]
datanew <- data.frame(data$SP, data$DP, data$SL, data$DL)
matrix <- datanew * vector
variance <- data.frame(var(matrix)[1], var(matrix)[2], var(matrix)[3], var(matrix)[4])


 var.matrix..1. var.matrix..2. var.matrix..3. var.matrix..4.
1    0.002225362   1.712515e-05   -0.001885247   -0.001308941

Thank you for your time, lmk if you want me to clarify anything.

CodePudding user response：

I can show you how to do what Nakagawa and Schielzeth describe in their paper, but I'm still not sure how it fits in with what you want to do.

Simulate data and fit example model:

library(lme4)
set.seed(101)
n <- 100
dd <- data.frame(x1 = rnorm(100), x2 = rnorm(100), x3 = rnorm(100),
                 f = factor(rep(1:10, each = 10)))
dd$y <- simulate(~ x1   x2   x3   (1|f),
       family = gaussian,
       newdata = dd,
       newparams = list(beta = c(1, 1, 2, 3), theta = 1, sigma = 1))[[1]]
m <- lmer(y ~ x1   x2   x3   (1|f), data = dd)

Here's what they say:

Estimating sigma^2_f can, in principle, be carried out by predicting fitted values based on the fixed effects alone (equivalent to multiplying the design matrix of the fixed effects with the vector of fixed effect estimates) followed by calculating the variance of these fitted values (Snijders & Bosker 1999).

However:

sigma^2_f is a single value, not a set of variances corresponding to variances explained by each fixed effect.
"multiplying" here refers to matrix multiplication (technically, matrix-by-vector multiplication)

## multiply fixed-effect model matrix X by FE coefficient vector
pred <- drop(getME(m, "X") %*% fixef(m))
## or equivalently:
pred2 <- predict(m, re.form = NA)
all.equal(pred, pred2) ## TRUE
var(pred)

(the last value might need to be multiplied by (n-1)/n to conform to S&N's "Note that sigma^2_f should be estimated without degrees-of-freedom correction.")

Note that R^2, variance proportions corresponding to fixed effects, and other related metrics for mixed models are a pretty deep statistical rabbit hole. I have used the r2glmm package to compute partial R^2 values for individual fixed effects within LMMs (not extended to GLMMs).

CodePudding user response：

It may interest you to check out the partR2 package, which was developed in part by the same people for getting fixed effect partial R2 values. You can read more about it from this paper on the package.

Here I give a simple demonstration with the carrots dataset included with lmerTest. First you load the requisite libraries: lmerTest for the mixed model and partR2 for the fixed effect R2 values. Then for simplicity's sake I dropped NA values and saved the data, as partR2 will have issues if you run the data as is since lmer will kick out missing values. Normally its wiser to handle missing data in a different way, but I'll skip that for now. Then you fit the model as you normally would. After, save the partR2 as an object. Here I've used 100 bootstraps but normally its better to use 1000 (its very slow, hence why I use 100 here). Finally I summarize the results at the end:

#### Load Libraries ####
library(lmerTest)
library(partR2)

#### Drop NA Values ####
carrots <- carrots %>% 
  drop_na()

#### Fit Model ####
fit <- lmer(Preference 
            ~ sens2 
              Homesize 
              (1| Consumer),
            data=carrots)

#### Extract Part R2 ####
part.fit <- partR2(fit,
       partvars = c("sens2","Homesize"),
       nboot = 100) 

#### Summary ####
summary(part.fit)

The summary can be seen here:

R2 (marginal) and 95% CI for the full model: 
 R2     CI_lower CI_upper ndf
 0.0643 0.04     0.0897   3  

----------

Part (semi-partial) R2:
 Predictor(s)   R2     CI_lower CI_upper ndf
 Model          0.0643 0.0400   0.0897   3  
 sens2          0.0498 0.0254   0.0754   2  
 Homesize       0.0145 0.0000   0.0425   2  
 sens2 Homesize 0.0643 0.0400   0.0897   1  

----------

Inclusive R2 (SC^2 * R2):
 Predictor IR2    CI_lower CI_upper
 sens2     0.0498 0.0317   0.0721  
 Homesize3 0.0145 0.0018   0.0315  

----------

Structure coefficients r(Yhat,x):
 Predictor SC      CI_lower CI_upper
 sens2      0.8804  0.7369   0.9838 
 Homesize3 -0.4749 -0.6765  -0.1797 

----------

Beta weights (standardised estimates)
 Predictor BW      CI_lower CI_upper
 sens2      0.2235  0.1782   0.2690 
 Homesize3 -0.2800 -0.4284  -0.0976 

----------

Side note: if you have a fixed effect included as a random slope, I believe partR2 cannot derive a fixed effect effect size. Don't remember why, but something to consider. Another thing to remember is that this package is fairly new (last year in fact), so what can be considered a "weak" or "strong" effect size is yet to be determined.