# Model assumption of linearity

I am trying to interpret the outcome of a test for assumption of linearity. This is the dataframe:

df <- structure(list(treatment = structure(c(1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("CCF",
"UN"), class = "factor"), random1 = structure(c(3L, 1L, 2L, 2L,
2L, 2L, 2L, 4L, 4L, 4L, 4L, 3L, 4L, 3L, 3L, 3L, 1L, 1L, 1L, 1L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), .Label = c("1.6",
"2", "3.2", "5", NA), class = "factor"), random2 = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L), .Label = c("1", "2", "3", "4", "5", "NA"), class = "factor"),
continuous = c(13.4834098739318, 13.5656429433723, 12.4635727711507,
18.72345150947, 18.4616104687818, 20.5685002028439, 13.8419704601596,
16.1418346212744, 17.2712407613484, 15.6206999481025, 17.3198253734436,
15.9326515550379, 13.6664227787624, 18.4006445221394, 15.9590212502841,
18.8509698995243, 20.5492911251772, 12.0971869009945, 14.2687663092537,
17.5558622926168, 12.0655307162184, 20.0060355952652, 15.9836412635937,
18.5999712367426, 14.9125382681618, 18.4091462029293, 18.766029822543,
15.8768079929326, 14.5894782578156, 11.6426318894049, 16.8206949611527,
17.0666712246649, 16.7071675430987, 16.2745705651548, 15.9203707655043
)), class = "data.frame", row.names = c(3L, 6L, 9L, 12L,
15L, 18L, 21L, 24L, 27L, 30L, 33L, 36L, 39L, 42L, 45L, 48L, 51L,
54L, 57L, 60L, 63L, 66L, 69L, 72L, 75L, 78L, 81L, 84L, 87L, 90L,
93L, 96L, 99L, 102L, 105L))


This in the model:

library(lme4)
model <- lmer((continuous) ~  treatment +(1|random1) + (1|random2), data= df, REML = TRUE)


Upon checking the relationship between the independent and dependent variables to be linear with Linearity<-plot(resid(model),df\$continuous) I got this result:

I was expecting my results to be either randomly scattered (linear relationship) or to show other behaviors (e.g. curvilinear relationship). How should I interpret this outcome almost completely on a straight line, and does it meet the linearity assumption?

This is what I get when I use Linearity<-plot(resid(model), fitted(model)):

This is the boxplot by treatments (colours are random1)

emmeans model results displaying estimated marginal means +- SE:

Example of random1 eeffect on CCF (different variable, see comments):

• It sounds like what you really wanted to look at was the plot of residuals versus fitted values, not the plot of residuals versus the dependent variable. In other words, plot(residuals(model) ~ fitted(model)). – Isabella Ghement Apr 8 '19 at 14:59
• thank you. Am I misinterpreting the the assumption that I am following from University of Illinois at Chicago or is it a typo? If I use your method I get something like the scatterplot above – BAlpine Apr 8 '19 at 15:49
• Also, random2 has only 2 levels, so you really do not want to fit random intercepts for this. – Robert Long Apr 8 '19 at 16:17
• And...... random1 has almost half it's values missing, corresponding to ALL of one of the two levels of random2....... – Robert Long Apr 8 '19 at 16:19
• Finally, please provide some detail about how these data were collected including the study design, and what the variables represent – Robert Long Apr 8 '19 at 16:20

Note that the model has a singular fit, and you should therefore not do any interpretation, or checking of assumptions, until that is resolved.

> model <- lmer(continuous ~  treatment + (1|random1) + (1|random2), data= df, REML = TRUE)

> isSingular(model)
[1] TRUE


There are some issues with your dataset:

1) Plot the data:

Straight away, we can see that there may be little hope of finding a significant treatment effect.

2) The factor random2 has only 2 levels, and this is not sufficient to warrant fitting random intercepts. The software will try to estimate a variance for a normally distributed variable with only 2 observations. This can not result in a meaningful estimate and you should not model this factor as random using lme4.

3) The factor random1 has 15 of 35 observations missing. Moreover these 15 observations correspond to ALL of the observations of a single random2 level.

Even after removing random2 as a random effect, there is still a singular fit. This is because there is virtually no intra-class correlation, so there is no need to fit random effects at all, with these data.

One way forward is to include random2 as a fixed effect - as if it were a potential confounder (so you would not try to interpret it's coefficient):

model1.6 <- lm(continuous ~  treatment + random2, data= df)

Residuals:
Min      1Q  Median      3Q     Max
-4.9496 -1.9370 -0.0953  1.9567  4.3314

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  16.2371     0.5454  29.771   <2e-16 ***
treatmentUN  -2.2660     1.8089  -1.253    0.219
random22      2.6211     1.8526   1.415    0.167

Residual standard error: 2.439 on 32 degrees of freedom
Multiple R-squared:  0.05887,   Adjusted R-squared:  4.796e-05
F-statistic: 1.001 on 2 and 32 DF,  p-value: 0.3788


It is hardly surprising that treatment is not significant, but this does handle the possible non-independence of observations due to clustering in random2. Note that we cannot do the same for random1 due to the extent of missing values.

• Thank you for your constructive comment. The design is unbalanced: 15 independent individuals for treatment UN and 20 for treatment CCF. random1 are areas within the CCF treatment. This is because only CCF was subject to an earlier experiment which I am not intereseted in , but I would like to account for in my model (it may have introduced variability). random2 was the batch of the instrument that I utilized to get my results.... – BAlpine Apr 10 '19 at 12:51
• ... What I would like to have is the difference between treatment. I am not sure about the point you discuss in random2: this is a random factor but each of my treatment have 15/20 individuals. I am not interested in random1 and random2, I just want to account for the variability that they introduce. Does it make sense? Also, I would appreciate if you can suggest a way to deal with the singular fit . Many thanks. – BAlpine Apr 10 '19 at 12:52
• @BAlpine I have updated the answer. You are not sure about my point about random2 ? What are you not sure about ? This factor has only 2 levels. If you fit random intercepts for it, the software will try to estimate a variance for a normally distributed variable that has only 2 observations - this is completely meaningless. – Robert Long Apr 10 '19 at 13:56
• Haha, OK. Well, that glmer model is also singular !! I don't think there is much more I can say about this. Comments are not intended for protracted discussion of substantive issues, and I have already edited the answer several times, following your comments. Please consider marking this as the accepted answer and if you have any further questions you can ask a new question.and I will be happy to to take a look at it. – Robert Long Apr 10 '19 at 15:41
• @BAlpine The last thing I would say on this is please do NOT use a singular model to produce plots or other outputs. if you would like to ask a new question please do so, but I don't have anything else to say on this one. These long discussions in comments with edits to the questions and answers can make the site difficult to use for future visitors. – Robert Long Apr 10 '19 at 16:02