I'm analyzing a three-way mixed linear model using lmer:

Y ~ Factor1 * Factor2 * Factor3 + (1|sensor)

However, different sensors have different gains, uniformly scaling the response of each sensor across conditions. Therefore, a multiplicative random effect seems more appropriate than the additive one I currently use. Can something like this be implemented in R?

Y ~ [Factor1 * Factor2 * Factor3] $\cdot$ (1|sensor)

(The dot stands for multiplication)

My dependent measure can be negative, so using log(Y) doesn't seem like a good solution.


I'd like to try a simplified, generalized formulation of this problem, using scalar notation. We'll start with a simple (fixed effects) linear model: $$y_i=X_i^1\beta_1+X_i^2 \beta_2+X_i^3 \beta_3+\epsilon$$ In this model, given $X$ and $y$, the $\beta$ coefficients can be easily estimated by ordinary least squares. Now, we'll complicate things by assuming that different subsets of observations are sampled by different sensors and each sensor has some random additive contribution $u$. This contribution is normally distributed with 0 expectancy and an unknown variance across sensors. $$y_i=X_i^1\beta_1+X_i^2 \beta_2+X_i^3 \beta_3+u_{sensor(i)}+\epsilon$$ This model can be estimated as a mixed effects model: Y~X1+X2+X3+(1|sensor). However, I'd like to consider the case of a multiplicative random effect: $$y_i=X_i^1\beta_1u_{sensor(i)}+X_i^2 \beta_2u_{sensor(i)}+X_i^3 \beta_3u_{sensor(i)}+..+\epsilon$$

Note that for each sensor, there's a single scalar gain. I wish to estimate this model by R. To my understanding, the standard random slope approach won't do it, since Y~X1+X2+x3+(X1+X2+X3|sensor) would estimate the model $$ y_i=X_i^1\beta_1u_{sensor(i)}^1+X_i^2 \beta_2u_{sensor(i)}^2+X_i^3 \beta_3u_{sensor(i)}^3+..+\epsilon $$ , estimating three different random gains for each sensor instead of one.Y~0+(X1+X2+X3|sensor) won't help either.

Any ideas?

  • $\begingroup$ I'm not sure I understand the question. Could you want Y ~ Factor1 * Factor2 * Factor3 + (Factor1 * Factor2 * Factor3|sensor)? $\endgroup$ – Roland Jun 4 '15 at 13:14
  • $\begingroup$ Please let me know if I'm wrong, the formula you stated will result with a model with a separate random effect for each condition (e.g. for a 3x7x2 design that would be 42 random effects). I'm interested in estimating only a single random effect, but a one that has multiplicative effect on y - instead of $$y=X\beta+u+\epsilon$$, I'd like $$y=X\beta\cdot u+\epsilon$$ $\endgroup$ – Trisoloriansunscreen Jun 4 '15 at 14:25
  • $\begingroup$ That is over my head, but it looks a bit like Y ~ 0 + (Factor1 * Factor2 * Factor3|sensor). $\endgroup$ – Roland Jun 4 '15 at 14:48
  • $\begingroup$ @Roland, I'm afraid that the similarity is only superficial. Please have a look at my edit above.. $\endgroup$ – Trisoloriansunscreen Jun 5 '15 at 10:10

There is an identifiability issue. Suppose $\beta$ and $u$ work. Then, $\beta/2$ and $2u$ will work as well.

The basic solution is to build standard curves for your sensors and run them with your experiments. Then, your raw measurements would be transformed into measurements on the same scale.

Failing that, you will need to let one sensor be the "reference sensor". The other sensors will have variability as a scale factor relative to the reference sensor.

Then, your model is nonlinear. Without linearizing it somehow you probably should not try to fit it using linear model methods. Instead, you should use some sort of nonlinear model method.

Here is an example using nonlinear least squares. In this example, I'm using continuous variates, though the same approach is useful for factors by using effects coding. It is just easier to demonstrate using continuous variates.

# Gin up some data.


Sensor <- data.frame(
    X1 = rnorm(150),
    X2 = rnorm(150),
    X3 = rnorm(150),
    U1 = c(rep(1, 50), rep(0, 50), rep(0, 50)),
    U2 = c(rep(0, 50), rep(1, 50), rep(0, 50)),
    U3 = c(rep(0, 50), rep(0, 50), rep(1, 50))

# Set some sensor gains and beta values.  Note
# that the first sensor is the reference sensor.

g <- c(1, 2, 4)
b <- c(-1, 1, 3)

# Simulate a simple version.  Careful of factors
# with numeric values.

Y <- with(
    (b[1]*X1 + b[2]*X2 + b[3]*X3) * (g[1]*U1 + g[2]*U2 + g[3]*U3)

e <- rnorm(150)

Sensor$Y <- Y + e

This duplicates the model you show above, with an additive error after the nonlinear transformation. Note this also assumes a zero intercept.

In this example, the first sensor is the reference sensor. The second and third sensors produce measurements that differ by a scale factor of 2 and 4 respectively from the measurement that the first sensor would produce.

# Fit a model using nonlinear least squares.


f <- function(b1, b2, b3, X1, X2, X3, g2, g3, U1, U2, U3) {

    (b1*X1 + b2*X2 + b3*X3) * (1*U1 + g2*U2 + g3*U3)


fit <- nls(
    formula = Y ~ f(b1, b2, b3, X1, X2, X3, g2, g3, U1, U2, U3),
    start= c(b1=1, b2=2.4, b3=7, g2=2, g3=3),
    data = Sensor


#Formula: Y ~ f(b1, b2, b3, X1, X2, X3, g2, g3, U1, U2, U3)
#   Estimate Std. Error t value Pr(>|t|)
#b1 -1.02810    0.06005  -17.12   <2e-16 ***
#b2  1.01436    0.06237   16.26   <2e-16 ***
#b3  3.12400    0.15062   20.74   <2e-16 ***
#g2  1.85447    0.09965   18.61   <2e-16 ***
#g3  3.86931    0.19169   20.18   <2e-16 ***
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Residual standard error: 1.106 on 145 degrees of freedom
#Number of iterations to convergence: 5
#Achieved convergence tolerance: 1.782e-07

Pretty good agreement all around with the $\beta$ coefficients, the sensor scale factors, and the error variance.


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.