# How to model this mixed model

I have following data:

Each column represents a specific manufacturer (that's why I skipped the title) and from each one I have several sensors. Now I would like to perform a mixed model. I thought about using the manufacturer and the sensors as random effects but I'm not quite sure about that.

In r resp. lme4 the syntax would be:

lmer( R ~ c + (Manufacturer | sensor), data = data)


Is this correct?

What you can see here is the response of the resistance of the sensors against several gas concentrations.

edit: More information: There are three manufacturers with two of 5 sensors each and one of 10 sensors. Basically it is:

Manufacturer A: SensorA1, SensorA2, SensorA3, SensorA4, SensorA5
Manufacturer B: SensorB1, SensorB2, SensorB3, SensorB4, SensorB5
Manufacturer C: SensorC1, SensorC2, SensorC3, SensorC4, SensorC5, SensorC6, SensorC7, SensorC8, SensorC9, SensorC10


All the sensors are housed in the same experimental unit and measured in parallel under varying conditions. Here in this case it is the applied concentration varying from 50 to 500 ppm. So all I did in the plot is to group these sensors according to their manufacturer. The colours represent the different concentrations (it is hardly visible in the index, I have to admit..), so basically the colours are not necessary as the axis already provides this information.

My goal now is to figure out how the sensors react to the different gases individually and how do the manufacturers.

Plots:

The formula

lmer(R~ poly(c, 1) + Manufacturer + (c | Manufacturer/Sensor), data = data)


results in

respectively in

and the formula

lmer(R~ poly(c, 1) + Manufacturer:c + (c | Sensor), data = data)


results in

I think the blue fit for the third group is not bad while the other two need an improvement. I'm playing a bit around but I'm not really able to enhance the situation. What I wonder is: Without any knowledge I would draw an almost straight line through the data points, separately for each Manufacturer (color). Why are above formulas not able to do so?

editedit: Matteo's suggestion lmer( I(log(R)) ~ I(log(c)) + Manufacturer + I(log(c)):Manufacturer + ( I(log(c)) | sensor), data = data) results in

I can't get my head around it why the fit behaves as it does. Why are they so off all the time?

my code:

lmer_model<- lmer(log10(R) ~ log10(c) + Manufacturer + log10(c):Manufacturer +
(log10(c) | Sensor), data = data)
ggplot(data, aes(x = c, y = R, colour = Manufacturer)) +
geom_point(alpha = 0.5) +
scale_x_continuous(
trans='log10'
, breaks = round(seq(min(data$$c), max(data$$c), by = 150),500)) +
scale_y_continuous(
trans='log10'
, labels=format_si()
, breaks = round(seq(min(data$$R), max(data$$R), by = 500000), 1000000)) +
theme_classic() +
geom_line(data = cbind(data, pred = predict(lmer_model)), aes(y = pred), size = 1) +
ylab(TeX("R / $$\\Omega$$ ")) + xlab(TeX("c / ppm")) +
theme(
legend.position = "none",
panel.spacing = unit(2, "lines"),
legend.text=element_text(size=7.5),
legend.box.margin=margin(-13,-13,-13,-13)
)

• I don't think that formula does what you want.. but just to be sure: what does sensor represent exactly? each dot is a measurement from a different sensor? (i.e., each of concentration levels is measured multiple times with different sensors?). Does the different manufacturers produce the "same" type of sensors? Dec 19, 2019 at 15:24
• I am having trouble understanding the experiment. Can you edit the question to give a bit more detail ? In particular is sensor nested in manufacturer ? Dec 19, 2019 at 22:26
• Firstly, thanks for the fast responses though I was absent for some time now. @matteo 'sensor' holds different sensors where each sensor represents all measurement points belonging to one specific sensor. The colours represent the different concentrations. In the plot you can only see three different manufacturers, from which there are 5 sensors from A, 5 sensors from B and 10 sensors from C.
– Ben
Jan 8, 2020 at 8:44
• @RobertLong Yes, sensor is nested in manufacturer. I will extend the background information!
– Ben
Jan 8, 2020 at 8:45

I think the model you want in lme4 syntax should be written as

lmer( R ~ c + Manufacturer + c:Manufacturer + ( c | sensor), data = data)


This model will have sensor-specific random intercept and slopes, and manufacturer specific (fixed) intercept and slopes, so I think it address your goal of estimating how individual sensors and manufactures respond to gas concentration.

I think that conceptually also Manufacturer should be a random effects, but because you have only three of them it is not very practical and not much to gain in trying to estimate the between-manufacturers variance.

One final note, by looking at your plot it seems that the relationship is linear when both variables (concentration and response) are expressed in logarithmic units. So I think it would be appropriate to log-transform the variables before entering them in the model. Note that this changes the interpretation of the coefficients which would code for relative difference instead of absolute differences as in standard linear regression. In a nutshell, a slope coefficient of $$\beta$$ would indicate that a $$x\%$$ increase in the predictor correspond approximately to $$x\beta\,\%$$ increase in the response, or more precisely to a change in the response equivalent to multiplying it by $$e^{\beta\log\left(\frac{100+x}{100}\right) }$$.

EDIT

I had a look at the data you posted, and fit the model I suggested above. The model seems to fits quite well, so there must be an issue in the way you computed or plotted the model predictions. I paste below the code and the resulting plot.

d <- read.table("data.csv",header=T,sep="\t")
library(ggplot2)
library(lme4)

# log transform variables
d$$log_c <- log(d$$concentration)
d$$log_R <- log(d$$Resistance)

# fit mixed model
m0 <- lmer( log_R ~ log_c + Manufacturer + log_c:Manufacturer + ( log_c | Sensor), data = d)

# not converging, try some more iterations
ss <- getME(m0, c("theta","fixef"))
m0 <- update(m0, start=ss, control=lmerControl(optCtrl=list(maxfun=2e4))) # OK

# make a new data frame for calculating model predictions
# (note that I am calculating them using only the fixed effects)
nd <- expand.grid(Manufacturer = unique(d$$Manufacturer), log_c = seq(min(d$$log_c),max(d$$log_c), length.out=100)) nd$$log_R <- predict(m0, newdata=nd, re.form=NA)
nd$$concentration <- exp(nd$$log_c)
nd$$Resistance <- exp(nd$$log_R)

# plotting
ggplot(data=d,aes(x=concentration, y=Resistance, color=Manufacturer))+geom_point()+ scale_x_continuous(trans = 'log10') +scale_y_continuous(trans = 'log10')+facet_grid(.~Manufacturer)+geom_line(data=nd)


• Thanks a lot for the great and helpful answer! I just noticed that you were aware of the log-transformations as well. I will take this into account now. Despite from that I would like to know why the number of Manufacturers plays a role in regards to the between-variances?
– Ben
Jan 8, 2020 at 14:21
• The point is that if you tell lme4 that Manufacturer is a random factor, e.g. with random manufacturer-specific intercepts, it will estimate the variance of these random intercepts. With only 3 manufacturer this would amount to trying to estimate a variance with only 3 data points: there is just not enough data to give a sensible estimate! Jan 8, 2020 at 15:20
• By the way, I think the discrepancy you see in the plots for the first two groups might be due to fitting the model on the original scale and then plotting in log-scale: try first transforming both R and c and then fitting the model. Also, the use of poly in this context seems superfluous (since you are asking to generate a polynomial of order 1, I think the only thing that the function does is centering and rescaling the predictor). Jan 8, 2020 at 15:25
• It doesn't matter but if you put log()in the lmer formula make sure to use the I() function, e.g. lmer( I(log(R)) ~ I(log(c)) + Manufacturer + I(log(c)):Manufacturer + ( I(log(c)) | sensor), data = data) Jan 8, 2020 at 19:31
• If you add them as interactions make sure you add them also as main effects. If you have many more possible predictions the number of interactions grows quickly and you might start overfitting, so for each terms that you add in the model you may want to check that it does not decrease the predictive ability of the model, using e.g. likelihood ratio test or cross-validation. Jan 13, 2020 at 8:06

In r resp. lme4 the syntax would be:

lmer( R ~ c + (Manufacturer | sensor), data = data)


Is this correct?

No. The above model specifies that observations are clustered within sensor (which is correct) but it also says that the effect of Manufacturer varies within each level of sensor (ie there are random slopes for Manufacturer), which is incorrect, since sensor is nested within Manufacturer. Also, you do not specify Manufacturer as a fixed effect which means that the random slopes vary around zero.

lmer( R ~ c + (1 | Manufacturer/sensor), data = data)


which is equivalent to

lmer( R ~ c + (1 | Manufacturer) + (1 | Manufacturer:sensor), data = data)


and assuming that sensor is coded uniquely as it seems to be in the information you added, it is also equivalent to:

lmer( R ~ c + (1 | Manufacturer) + (1 | sensor), data = data)


These models fit fixed effects for concentration and random intercepts for Manufacturer and sensor. The main issue I see with this is that there are only 3 levels of Manufacturer so you are unlikely to get a good estimate of it's variance. So an alternative is to specify Manufacturer as a fixed effect:

lmer( R ~ c + Manufacturer + (1 | sensor), data = data)


You might also want to allow the effect of concentration to vary across sensors by specifying random slopes for concentration:

lmer( R ~ c + Manufacturer + (c | sensor), data = data)

• Thank you very much for the detailed line-to-line explanations. This helps me a lot to dive deeper into the topic! Here I have, among others but as I would like to learn the concepts, to get started the first question why the levels of a variable play a role when it comes to treat them as random effects?
– Ben
Jan 8, 2020 at 14:23
• You're welcome. Do you mean the number of levels ? Jan 8, 2020 at 15:21
• That's because the random effects are assumed to be normally distributed and the software will estimate a variance for each variable, so in the case of your manufacturers you only have 3 levels so it means the software has to estimate the variance from only 3 observations. You can still do it, but the estimate is unlikely to be very good and you have no way of knowing how good it is. Jan 8, 2020 at 15:42
• Yes. Imagine that you wanted to estimate the variance of exam scores for 200 student's taking the same exam. If you only took 3 samples then you are unlikely to get a good estimate. Jan 9, 2020 at 15:05
• Kind of but I can't manage to apply the model description successfully, hence, I would prefer to keep it open.
– Ben
Jan 10, 2020 at 16:34