2
$\begingroup$

I have a question regarding mixed linear effects model.
For example, now I have 100 patients, and each patient has generated 10 to 30 eggs. So the dataset has 2500 rows (one egg per row). The response variable 'mature' is a binary one, 1 means this egg matures and 0 means not mature. We are interested in a variable called 'genotype' which is a trait per patient. So all the rows of same patient will have same genotype, but different patient could have same genotypes. There are 9 different genotypes. We also have some covariates, age and BMI.
So I want to see if genotypes associate with maturation of eggs. I know several options:

  1. Logistic regression: This one works but it assumes all eggs are independent, but we know eggs from same patient relate to each other.
  2. Generalized Mixed effects model: Similar to 1 but add a random effects of patient. I am confused about this one because I worried that the random effects(patient) would take away variation which could be explained by genotype, because genotype is the same to certain patient. I am not sure about this.


So this is my question: can I use Generalized Mixed effects model in this case to check relationship between genotype and egg maturation? If I can, which means the random effects is not designed to take away variance which should go to any fixed effects, how should I test the coefficient on genotype in this case, anova? If not, what is the best solution to this question? Aggregate the data to one patient one row and do Generalized linear model with beta distributed response (maturation rate), and put number of eggs per patient as weight?

Thank you for any input.

$\endgroup$
1
$\begingroup$

First, a comment: a logistic regression is totally fine for both your cases 1 and 2, and I think the question is whether to include person-level intercepts, and whether they are random or fixed effects, rather than what distribution the response should have. You've got a collection of binary random variables so you can't be wrong to say each one is Bernoulli distributed.


Let $y_{ij}$ be the measurement for egg $j$ from person $i$, where $i = 1, \dots, m=1000$ and $j = 1, \dots, n_i \leq 30$, and let $x_i \in \mathbb R^{11}$ give the person-level covariates for the 9 genotypes along with age and BMI.

It seems that you're comparing the following two models:

$\newcommand{\ilogit}{\text{logit}^{-1}}$ $$ P(y_{ij} = 1) = \ilogit(\mu + \alpha_i + x_i^T\beta) \tag{M1}\label{1} $$ where $\alpha \sim \mathcal N(0, \sigma^2_\alpha I_m)$ is a random intercept, and $$ P(y_{ij} = 1) = \ilogit(\mu + x_i^T\beta) \tag{M2}\label{2} $$ where now there is no person-level effect at all. Also really these are both conditioned on $x_i$ but I'm omitting that for simplicity.

$\ref{1}$ allows for correlation between eggs for a particular person, while $\ref{2}$ ignores this and says that when $x = 0$ we should predict the same thing for everybody, rather than having a unique baseline prediction per person.

Without even going to the statistics of the model, is $\ref{2}$ really reasonable? This should be your main concern. We know in general that there's variation in stuff like this between people, so scientifically I think $\ref{2}$ is probably a bad model. You want to model the things you know, and if one of those things is that the eggs from a given person are likely to be correlated, then you should explicitly account for this.

Now statistically, here's a heuristic way to picture what adding $\alpha$ does to your model. Our linear predictor is $\eta = \mu\mathbf 1 + X\beta + Z\alpha$ and $$ Z = \mathbf 1_{n_1} \oplus \dots \oplus \mathbf 1_{n_m} =\left(\begin{array}{ccccc}1 & 0& 0 & \dots & 0 \\ 1 & 0 & 0& \dots & 0 \\ \vdots & \vdots & \vdots & \dots & \vdots \\ 0 & 1 & 0 & \dots & 0 \\ 0 & 1 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \dots & \vdots \\ 0 & 0 & 0 & \dots & 1 \\ 0 & 0 & 0 & \dots & 1 \end{array}\right). $$

If we were doing a fixed effects linear regression, we could picture this as follows: in order to estimate $\beta$, first we project $y$ into the space orthogonal to $Z$. This means we're effectively centering each person's egg measurements. Then for each coefficient $\beta_j$ we look at what is explained by that variable after we've further made $y$ orthogonal to the other 10 predictors.

With our random intercept the idea is similar, although because of the shrinkage we no longer exactly have orthogonal projections. But I think it's fair to still think of $Z\alpha$ as effectively centering each person, and i think that gives us a better measurement of the fixed effect trends, because otherwise we might be thinking there's a trend when it's really just caused by a (partial) correlation between age, say, and the baseline.

Maybe a succinct way to answer the question in your title is that adding a predictor will increase the variance explained, but a fixed effect like age attempts to describe the part of the data orthogonal to the other predictors, so we sometimes get a much clearer picture by deliberately removing variation.

As to whether the $\alpha_i$ should be treated as random or fixed, I don't think anyone would argue that they make sense as a fixed effect, so I'm not going to discuss that further.


Update

In response to your comments: regarding your question 1, I go through the math of what i'm talking about here. By adding an additional predictor you probably will indeed explain more of total variation and $R^2$ will probably increase. But that isn't necessarily good or bad. What matters for a single predictor is the explained variation after the other predictors are accounted for, and that's what I show in that linked answer of mine. So adding $Z$ we are now forcing age to only explain the part of the variation that is not explained by the other predictors, which now includes a person-level intercept. But in cases like this that can be a good thing.

Here's an example: we've got a linear regression with two predictors, a continuous x and a binary grp. In this case it turns out that x is correlated with grp in that when grp=2, x is larger. But y is actually only a function of grp.

Here's what the data look like. You can see that within a group there's no relationship between x and y, but due to the correlation between x and grp if we just do the model y ~ x then the slope will be significant, i.e. we will think we've learned something. But really x is confounded with grp, and once we've accounted for grp, as we do when we look at the residuals of the model y ~ grp, we see then that x explains basically nothing of the remaining variation.

ex

That's why it's so important to include person-level intercepts if you think there's a reason for them to be there: if you don't include them you may be inviting confounding between your other predictors and the group level effects which are still there, just unmodeled. This can lead you to have confounded effects appear significant which would be a pretty serious mistake if it could have so easily been avoided.

Here's the code for the example:

set.seed(123)
dat <- data.frame(
  grp=rep(1:2, each=25),
  x=rnorm(50, rep(c(0,10), each=25))
)
dat$y <- dat$grp + rnorm(50, 0, .25)


summary(m.full <- lm(y~., data=dat))
summary(m.x <- lm(y~x, data=dat))

par(mfcol=c(1,2))
plot(y~x, data=dat, col=ifelse(dat$grp==1, "red", "blue"), pch=19, main="Original Data")
legend("topleft", c("group 1", "group 2"), col=c("red", "blue"), pch=19, bty="n")
legend("bottomright", "regression line for y~x", bty='n', lwd=2)
abline(m.x, lwd=2)

m.grp <- lm(y~grp, data=dat)
resid.grp <- resid(m.grp)

plot(resid.grp ~ x, data=dat,col=ifelse(dat$grp==1, "red", "blue"), pch=19, main="Residuals of y ~ grp versus x")
legend("bottom", c("group 1", "group 2"), col=c("red", "blue"), pch=19, bty="n")
legend("top", "resid ~ x", bty="n", lwd=2)
abline(lm(resid.grp ~ dat$x), lwd=2)

par(mfcol=c(1,1))
$\endgroup$
  • $\begingroup$ Thank you very much for your detailed answer! Two questions: 1. When you are saying 'beta is estimated after making y orthogonal to Z', if I understand this right, does it mean that the part of variance in y which could be explained by Z or age individually, now will be explained by Z? So are we underestimating age's effect in this case? This is what I worried. $\endgroup$ – Li Sun Apr 30 '18 at 21:40
  • $\begingroup$ Sorry the second question is, if use fixed effect for patient, does this use a lot more degrees of freedom and lower the power of this model? $\endgroup$ – Li Sun Apr 30 '18 at 21:46
  • $\begingroup$ @LiSun i've just added an update to better explain what i mean by that orthogonalizing $\endgroup$ – jld May 1 '18 at 3:10
  • $\begingroup$ @LiSun regarding your second comment, fixed effects do indeed use more degrees of freedom because they're not shrunk. I don't really work with power computations so i don't want to answer that without some more thought and say something wrong but i suspect you're right. I personally generally prefer random effects when appropriate because they have better predictive ability. The estimates are much more stable for groups with small numbers of observations and i often have factors with lot of levels, and you can make predictions for levels you never saw when training which is pretty neat $\endgroup$ – jld May 1 '18 at 3:16
  • $\begingroup$ Thank you very much for your patience and nice presentation of your answer. I totally agree with you about your opinion. However, what if, in your example, there is only 1 possible 'x' value in each 'grp'? Which is the case in my study, 1 patient('grp') could have only 1 genotype('x'). In this case, the 'real effects' of x could be attributed to grp, because they are perfectly correlated now. So my question is in this case, what should I do. $\endgroup$ – Li Sun May 3 '18 at 15:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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