# Linear mixed model beta coefficient larger than one

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

and here is the figure of the relationship between algae1 and protist1

• 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.
– EdM
Commented Jul 2, 2023 at 16:19
• sorry, the coefs function is form piecewiseSEM and the VIF function is from semEff package Commented Jul 2, 2023 at 16:49
• 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 ... Commented Jul 2, 2023 at 18:53
• 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)) Commented Jul 2, 2023 at 19:16
• @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? Commented Jul 6, 2023 at 7:27

• @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."
• 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 ... Commented Jul 2, 2023 at 17:08