2
$\begingroup$

My question is related to the linear mixed model. I have data on the built land area for three sampling sites from which we collected samples for 11 consecutive days. We collected the community of bacteria, algae, and protists using a DNA-based method. Our goal is to examine the relationships among all variables.

Since all of these data have high dimensionality, we used principal component analysis to reduce the dimensionality of each group of data. Subsequently, we used the first axis of each variable to represent the changes in bacteria, algae, and protists.

Based on ordination analyses, I observed that the data collected from each site were clustered closer together. Therefore, I believe that a mixed model, taking into account the sampling site as a random effect, would be appropriate. Additionally, the diagnostic plot for the independence of samples suggested a pseudo-replication problem arising from the same sampling site. For these reasons, I constructed several linear mixed models using the lmer function from the lme4 package.

However, while diagnosing each individual model, I noticed something that seemed a bit off. The code below illustrates my confusion. So, I built three models as follows:

library(lme4)
library(dplyr)
library(piecewiseSEM)

lme4::lmer(protist1 ~  pro1   +
      (1 | Site) , data = data1) %>% coefs

the estimate of pro1 was 0.4483 while the standardized estimate was 1.3535 ***

lme4::lmer(protist1 ~  algae1   +
         (1 | Site) , data = data1) %>% coefs

the estimate of algae1 was 0.65 while the standardized estimate was 0.899 ***

But when I run the following model:

lme4::lmer(protist1 ~  pro1   + algae1 + 
         (1 | Site) , data = data1) %>% coefs

the estimate of pro1 was 0.137 while the standardized estimate was 0.415; the estimate of algae1 was 0.48 and the standardized estimate was 0.668***

lme4::lmer(protist1 ~  pro1 + algae1   +
         (1 | Site) , data = data1)%>% VIF 

the VIF value for pro1 and algae1 was 2.79. which I assumed was okay.

I understood that the standardized estimate can sometimes be larger than one when a model includes multiple predictors, which could be attributed to high collinearity among the predictors, especially in a linear model. However, in this case, the situation was reversed. In the first example, the standardized estimate of pro1 was larger than one, even exceeding the non-standardized estimate, which has left me confused.

From what I understand, the standardized estimate is obtained from scaled data in R, and it is not expected to be larger than +1 or smaller than -1. I might be missing something in the concept of the mixed model, or perhaps using PCA components has an impact on the scaling function and subsequently on the estimates and standard estimates. How can I explain the findings in these models?

dataset

data1 <- data.frame(pro1 = c(13.384, 7.738, 16.957, 13.879, 4.038, 18.585, 15.204, 15.021, 12.217, 14.731, 6.054, -1.077, -0.669, 4.245, -1.841, 0.767, 11.68, 14.147, 11.694, 13.767, 15.845, -16.768, 
-18.923, -22.405, -20.059, -19.307, -17.848, -15.305, -17.387, -19.718, -19.168,  -19.477),
algae1 = c(9.73, 2.698, 3.08, 1.885, 3.677, 7.5, 6.923, 7.446, 7.141, 8.472, -2.232, -2.308, -2.178, -4.056, -6.327, -2.901, 7.754, 8.303, 7.983, 5.946, 9.527, -8.026, -8.31, -8.355, -9.927, 
-8.161, -6.754, -5.298, -4.841, -5.848, -4.499,  -8.045),
protist1 = c(5.201, 3.667, 2.76, 1.186, 0.976, 6.545, 2.279, 3.279, 2.935, 3.991, 0.104, -0.787, 
1.925, -2.069, -4.589, -2.215, 5.698, 7.084, 5.588, 7.445, 8.541, -4.416, -3.863, -6.164, -5.991, -6.039, -6.09, -4.103, -4.023, -5.758, -6.377, -6.719),
Site = c("W1", "W1", "W1", "W1", "W1", "W1", "W1", "W1", "W1", "W1", "W2", "W2", "W2", "W2", 
"W2", "W2", "W2", "W2", "W2", "W2", "W2", "W3", "W3", "W3", "W3", "W3", "W3", "W3",
                    "W3", "W3", "W3",  "W3"))

here is the figure of the relationship between pro1 and protist1 enter image description here

and here is the figure of the relationship between algae1 and protist1 enter image description here

$\endgroup$
5
  • $\begingroup$ Please edit the question to specify the package that provides coefs and VIF. Although the example data are very helpful, I can't reproduce all your results without that information about coefs and VIF. $\endgroup$
    – EdM
    Commented Jul 2, 2023 at 16:19
  • $\begingroup$ sorry, the coefs function is form piecewiseSEM and the VIF function is from semEff package $\endgroup$
    – Bob Adyari
    Commented Jul 2, 2023 at 16:49
  • 1
    $\begingroup$ I don't have enough for an answer yet. The short version of the answer is "the relationship (linear regression coefficient based on (y/sd(y)) ~ (x/sd(x)) == corr(y,x) doesn't hold for linear mixed models". I was going to mumble something about "the formula for the univariate regression coefficient, cov(y,x)/sd(x), reduces to cor(x,y) when sd(x) = sd(y) = 1 when we scale the predictors, but in this case what we're computing is something akin to lm(r1 ~ pro1) where r1 is the residuals from lm(protist1 ~ Site). That's not quite right though ... $\endgroup$
    – Ben Bolker
    Commented Jul 2, 2023 at 18:53
  • $\begingroup$ In this case (with only three sites) I might recommend a fixed effect of size, and allowing for site-by-covariate interactions: options(contrasts = c("contr.sum", "contr.poly")); coef(lm(protist1 ~ (pro1*Site), data1_sc)) $\endgroup$
    – Ben Bolker
    Commented Jul 2, 2023 at 19:16
  • $\begingroup$ @BenBolker, actually, there are more variables in my study and my goal was to construct a mixed-model-based SEM using piecewiseSEM package and I only want focus on the variation that can be explained by my variables of interest. However, I do not quite get how the model would look when the site was considered a fixed effect in linear model-based SEM. I read from some references that 3 levels of random effect might not be sufficient but with 11 days (replicate) would it still be not okay? and also, is there any concern regarding the coeffcient larger than one in single predictor mixed model? $\endgroup$
    – Bob Adyari
    Commented Jul 6, 2023 at 7:27

1 Answer 1

6
$\begingroup$

Welcome to the site. You wrote:

From what I understand, the standard estimate is obtained from scaled data in R, and it is not expected to be larger than one or smaller than negative one

AFAIK, this understanding is incorrect. The standard estimate is, indeed, obtained from scaled data, but it does not follow that it should be bounded by -1 and 1. It just means that the parameter estimate is now applied to the scaled data. So, if the standard estimate is (say) 5, that means that the DV goes up by 5 units for each 1 SD increase in the IV.

$\endgroup$
4
  • 1
    $\begingroup$ @BenBolker the OP used the coefs() function from the piecewiseSEM package to accept the pipe from lmer() . That extracts by default standardized coefficients such that "Raw coefficients are scaled by the ratio of the standard deviation of x divided by the standard deviation of y." $\endgroup$
    – EdM
    Commented Jul 2, 2023 at 16:54
  • 1
    $\begingroup$ @PeterFlom it was not SE, it was the standardized estimated that was calculated using the ‘coefs' function from piecewiseSEM package. So as said before, the standardized estimate was larger than one. I tried to understand this from linear model concept. Isn't the std estimate for one predictor same as the bivariate correlation coefficient? and if so, in linear model, it was impossible to get correlation larger than one. However, I do not quite understand when it is the case of linear mixed model. I guess it is not quite the same but what is the intuitive explanation for that? $\endgroup$
    – Bob Adyari
    Commented Jul 2, 2023 at 16:55
  • $\begingroup$ @BenBolker Should we just delete all the comments? $\endgroup$
    – Peter Flom
    Commented Jul 2, 2023 at 17:06
  • 2
    $\begingroup$ Probably. I now see the OP's point, the short answer is "the relationship (linear regression coefficient based on (y/sd(y)) ~ (x/sd(x)) == corr(y,x) doesn't hold for linear mixed models", but I'm trying to come up with an explanation now ... $\endgroup$
    – Ben Bolker
    Commented Jul 2, 2023 at 17:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.