2
$\begingroup$

I have a dataset containing one dependent variable which is the concentration of antibiotic needed to kill a bacteria, which was measured for several different antibiotics for three different microorganisms. The antibiotics are also divided in two groups based on their origin (synthetic or natural).

The data can be described as follows:

 $ ID: Factor w/ 3977 levels "1","2","3","4",..: 4 5 9 10 11 12 13 14 15 16 ...
 $ OR: Factor w/ 2 levels "natural", "synthetic": 2 2 2 2 2 2 2 2 1 2 ...
 $ MC: Factor w/ 3 levels "M1","M2","M3": 1 1 1 1 1 1 1 1 1 1 ...
 $ Y : num  1.745 0.125 2.301 -1.615 -2.026 ...

Additionally, as you can see the dataset is quite unbalanced and as a lot of missing values.

                    MR
OR          M1      M2        M3
natural   1267    1032       400
synthetic 2129    2044       944

I have specified a couple of formulas for the lmer() model.

(a) Y ~ OR * MC + (1|ID) 
(b) Y ~ OR + MC + (1|ID)
(c) Y ~ OR + MC + (OR+MC|ID)

For model (a), Anova with type 3 error showed that OR:MC is not significant.

Model (b), shows a slope on the residuals, so i tried model c.

Model (c) does not run in R (Error: number of observations (=7816) <= number of random effects (=27839)) so i turned to matlab (also runs on julia), and also shows the residuals to have a slope.

The slope in the residuals can be attributed, from what i understand, to several issues, poorly specified random effects or autocorrelation. The fact is that there might be autocorrelation as some antibiotics differ from other antibiotics in just a few atoms.

Any idea on how to properly specify the model?

Edit:

y=residuals; x=fitted residuals vs fitted

y=residuals; x=observed residuals vs observed

Model d (with a random slope and intercept for all levels of OR:MC)

(b) Y ~ OR * MC + (OR:MC|ID)

I believe both model b, c and now d are well specified, model b as a logLik of -7981, c of -7944 and d of -7933. Suggesting d is the better.

Improve this question
$\endgroup$
6
  • $\begingroup$ The problem probably lies in the error message from R so I would focus on this and not try to get around it by running the same model in Matlab. Which model gave this error ? $\endgroup$ – Robert Long Jul 1 '16 at 18:26
  • $\begingroup$ I'm betting that the error comes from the model with 2 random slopes. Basically you don't have enough observations to enable the estimation of the random slopes. Stick with the random intercepts model. Btw, what exactly is ID ? $\endgroup$ – Robert Long Jul 1 '16 at 18:51
  • $\begingroup$ I will edit the end part of the question as i think the a) b) and c) were misunderstood, i am sorry. These refer to the models. ID is the subject (antibiotic) identification. $\endgroup$ – SwatchPuppy Jul 2 '16 at 0:33
  • $\begingroup$ Random intercepts also generate models with random slopes on the residuals. $\endgroup$ – SwatchPuppy Jul 2 '16 at 0:36
  • $\begingroup$ OK, it's as I suspected, model c) isn't viable because you don't have enough observations (hence the error), What do you mean "shows a slope on the residuals" - can you post the residual plot ? What exactly is ID (what are the clusters ?) $\endgroup$ – Robert Long Jul 2 '16 at 5:13
3
$\begingroup$

Model c) is not viable because you do not have enough observations.

In model b) a linear relationship (slope) of the plot of residuals vs Y is to be expected because you have fitted a linear (mixed) model. Taking the linear mixed model formulation:

$$\mathbf{Y}=\mathbf{X \beta} + \textbf{Zb} + \epsilon,$$

where $\mathbf{X}$ is the model matrix for the fixed effects, $\mathbf{\beta}$ is the fixed effects coefficient vector, $\textbf{Z}$ is the model matrix for the random effects, $\textbf{b}$ is the random effects vector and $\epsilon$ is the error term vector, this can be trivially rearranged as:

$$ \epsilon = \mathbf{Y} - ( \mathbf{X \beta} + \textbf{Zb} ) $$

Since the residuals can be thought of as estimates of the errors, it follows that increasing values of $\mathbf{Y}$ , will be associated with increasing residuals. This can be seen in the following example:

> require(lme4)

> m0 <- lmer(Reaction ~ 1 + Days + (1+Days|Subject), sleepstudy)
> res <- resid(m0)
> plot(sleepstudy$Reaction, res)

enter image description here

Better diagnostic plots are a qq plot of the residuals and a plot of residuals vs fitted values:

> qqnorm(res)

enter image description here

> fits<- fitted(m0)
> plot(fits, res)

enter image description here

Improve this answer
$\endgroup$
5
  • $\begingroup$ I made a mistake, the plot i show if for the residuals vs fitted. And this is the one with the slope. $\endgroup$ – SwatchPuppy Jul 2 '16 at 10:17
  • $\begingroup$ Oh :( Well, the slope in the residuals vs fitted plot may indicate a missing variable - how does the plot look when you include the interaction as a fixed effect ? $\endgroup$ – Robert Long Jul 2 '16 at 12:06
  • $\begingroup$ Exactly the same. The only way for having the residuals without a slope, and i think it's wrong, is to define (ID|MC:OR) or any combination of (ID|MC) or (ID|OR). $\endgroup$ – SwatchPuppy Jul 2 '16 at 19:46
  • $\begingroup$ Please post the model output from Y ~ OR * MC + (1|ID) along with the qq plot of residuals. $\endgroup$ – Robert Long Jul 3 '16 at 10:52
  • $\begingroup$ The qq-plot seems to be normal. $\endgroup$ – SwatchPuppy Jul 4 '16 at 22:40

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.