How to model my data with linear mixed models for contrasts analysis 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

y=residuals;  x=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.
 A: 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)


Better diagnostic plots are a qq plot of the residuals and a plot of residuals vs fitted values:
> qqnorm(res)


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


