Test for linearity between two variables

The title might not be enough to understand my problem but i really did not know how to make it short but complete. My question might be silly but, it is also silly to stay stuck too long, so here I go: I have a set of growth and dilution values and I want to know if there is significant linear relationship between these two variables to be able to get the slope and the y-intercept and later compare them with other.

   dilution net_growth
1       1.0      1.263
2       1.0      0.047
3       1.0      0.088
4       0.6      1.692
5       0.6      0.447
6       0.6      0.671
7       0.3      2.271
8       0.3      1.307
9       0.3      1.485
10      0.1      3.338
11      0.1      2.333
12      0.1      2.141


Here is an example of the data and the plot resulting: plot(dilution, net_growth)

I wanted to simply use lm() to test the linear relationship and then an anova to compare the slope, but I got lost in all the different answers when looking for a solution on the net. With all the assumptions of the linear regression and the assumption of linearity itself, the qqplot that can help you test these assumptions but that can be interpreted in different ways. I think I mix everything.

So in short what can i do to test the linear relationship between my two variables and thus get a slope and y-intercept that i can use later.

• summary(lm(net_growth ~ dilution, yourData)) – PoGibas Oct 18 '17 at 9:48
• To add to @PoGibas comment, save the result of the regression fit <- lm(...) and use it later. summary gives you lots of information, but there's also, for instance, coef, residuals, etc. See the help page ?lmfor more. – Rui Barradas Oct 18 '17 at 9:52
• "I want to know if there is significant linear relationship between these two variables" For that you don't even need the lm model. You can simply do cor.test(DF$dilution, DF$net_growth). Note how that gives the same p-value as the linear model in Eric's answer. – Roland Oct 18 '17 at 11:48
• @Roland, it seems like Justine also would like to obtain "a slope and y-intercept that i can use later", e.g. compare slopes and intercepts for different cases. – Sextus Empiricus Oct 19 '17 at 8:07
• @MartijnWeterings I know. My comment is partly for the benefit of others who find this question through its title. – Roland Oct 19 '17 at 8:11

Why do you want to test for a significant linear fit?

Your data (as shown in your plot) passes the IOTT: The interocular trauma test. It hits you between the eyes. With such a small N (12 points) you might not get a very low p but ... so what?

In addition, you probably have substantive grounds for thinking that there is a linear relationship. If so, then you would mostly be looking for evidence against linearity in your own data - and that would be a sign of a problem in your data.

Looking at your data, I'd be more interested in trying to figure out why there is so much variation in growth at each level of dilution.

• And I would wonder why there seems to be a hint of the plotted values being bounded below at zero. – mdewey Oct 18 '17 at 12:37
• Looking at the plot, I'd bet the assumption of independence is violated and OP has repeated measures. If that is the case they should account for that dependency structure in their model and possibly use a linear mixed effects model. – Roland Oct 19 '17 at 6:34
• Exactly what conclusions "hit you between the eyes" in this plot, Peter? It's unclear even what you mean by "linear fit:" are you referring to fitting these data with a linear function or, perhaps, to assessing the possible departures from such a fit? – whuber Feb 6 '20 at 16:13

This answer was posted before the question was migrated to CV. A concise illustration of the assumptions in the linear regression model, illustrated with R code and data, can be found here.

df <- data.frame('dilution' = c(1, 1, 1, 0.6, 0.6, 0.6, 0.3, 0.3, 0.3, 0.1, 0.1, 0.1),
'net_growth' = c(1.263, 0.047, 0.088, 1.692, 0.447, 0.671, 2.271, 1.307,
1.485, 3.338, 2.333, 2.141))

fit <- lm(net_growth ~ dilution, df)

summary(fit)
#>
#> Call:
#> lm(formula = net_growth ~ dilution, data = df)
#>
#> Residuals:
#>     Min      1Q  Median      3Q     Max
#> -0.7478 -0.4281 -0.1947  0.4167  0.9992
#>
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)   2.5676     0.3176   8.083 1.08e-05 ***
#> dilution     -2.2880     0.5258  -4.352  0.00144 **
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 0.6176 on 10 degrees of freedom
#> Multiple R-squared:  0.6544, Adjusted R-squared:  0.6199
#> F-statistic: 18.94 on 1 and 10 DF,  p-value: 0.00144

plot(net_growth ~ dilution, df);  abline(lm(net_growth ~ dilution, df))


With ggplot2 you can overplay the 95% confidence interval to the plot, see here for more options,

# install.packages(c("ggplot2"), dependencies = TRUE)
library(ggplot2)
ggplot(df, aes(dilution, net_growth)) +
geom_point() +
geom_smooth(method=lm)


• Thanks a lot. So there is no need to check any assumption at this step ? – Justine Courboules Oct 18 '17 at 10:32
• @JustineCourboules The answer above is under the assumptions that underly this kind of linear model, see this link for more details. You can check you assumptions by calling plot on the lm object. This will yield a number of diagnostic plots. – Paul Hiemstra Oct 18 '17 at 10:37
• Your assumptions is crucial! In addition to the link provided by @PaulHiemstra I want to point to this handy tutorial on the linear regression model illustrated with R code and data. – Eric Fail Oct 18 '17 at 10:47
• @PaulHiemstra Wait I misundertsand your first answer actually. The p-value given in the summary of the lm() is accurate only under the assumptions of the linear model? But then is there a way to check them other than with the diagnostic plots? – Justine Courboules Oct 18 '17 at 11:03
• @JustineCourboules in some cases the check is very clear, e.g. in the case of using a normal distribution for clearly log-normally distributed data. When things get closer to the acceptable/non-acceptable line it is much harder to provide objective boundaries that can be applied without context. There you will have to make a value judgement yourself. – Paul Hiemstra Oct 18 '17 at 12:53

What I noticed is that there seems to be a correlation between your residual terms and the order of measurement (and this is also rather typical for many ordered experiments done by others, who often ignore it as well).

The first measurement is always higher than the other two. This indicates that there may be some random variation in the intercept, while the slope might be very precise each time (each order). (variations of the theme exist, such as variations in the slope with every change in the experiment order)

If you incorporate this order in the analysis then you may obtain a higher precision for the estimate of the slope.

See the below image and code for the difference. I have also added a quadratic term which appears to work well.

# data
inp <- c(1.0,1.263,1.0,0.047,1.0,0.088,0.6,1.692,0.6,0.447,0.6,0.671,0.3,2.271,0.3,1.307,0.3,1.485,0.1,3.338,0.1,2.333,0.1,2.141)
x <- inp[(1+c(0:11)*2)]
y <- inp[(2+c(0:11)*2)]
trial <- c(1,2,3,1,2,3,1,2,3,1,2,3)

# plot
plot(x,y,pch=21,bg=1,col=1)
lines(x[which(trial==1)],y[which(trial==1)],col="gray",lty=2)
lines(x[which(trial==2)],y[wfunction hich(trial==2)],col="gray",lty=2)
lines(x[which(trial==3)],y[which(trial==3)],col="gray",lty=2)

# model 1
model1 <- lm(y~1+x)
summary(model1)
lines(x,predict(model1),lwd=2,col=2)

# model 2
library(lme4)
model2 <- lmer(y ~ 1 + x + (1|trial))
summary(model2)
lines(x[which(trial==1)],predict(model2)[which(trial==1)],col=3)
lines(x[which(trial==2)],predict(model2)[which(trial==2)],col=3)
lines(x[which(trial==3)],predict(model2)[which(trial==3)],col=3)

# note that the an additional quadratic term is a statistical significant term
x2 <- x^2
model2b <- lmer(y ~ 1 + x + x2 + (1|trial))

anova(model2,model2b)

newdata <- as.data.frame(list(y=c(10:100)/100,x=c(10:100)/100,x2=(c(10:100)/100)^2,trial=rep(1,91)))
lines(newdata$x,predict(object=model2b,newdata=newdata),col=4) newdata <- as.data.frame(list(y=c(10:100)/100,x=c(10:100)/100,x2=(c(10:100)/100)^2,trial=rep(2,91))) lines(newdata$x,predict(object=model2b,newdata=newdata),col=4)
newdata <- as.data.frame(list(y=c(10:100)/100,x=c(10:100)/100,x2=(c(10:100)/100)^2,trial=rep(3,91)))
lines(newdata\$x,predict(object=model2b,newdata=newdata),col=4)


pieces from the output of the above:

model 1

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   2.5676     0.3176   8.083 1.08e-05 ***
x            -2.2880     0.5258  -4.352  0.00144 **


model 2 has lower standard error for the linear term (0.2552 instead of 0.5258)

Random effects:
Groups   Name        Variance Std.Dev.
trial    (Intercept) 0.36454  0.6038
Residual             0.08986  0.2998
Number of obs: 12, groups:  trial, 3

Fixed effects:
Estimate Std. Error t value
(Intercept)   2.5676     0.3812   6.736
x            -2.2880     0.2552  -8.966


> anova(model2,model2b)
Models:
model2: y ~ 1 + x + (1 | trial)
model2b: y ~ 1 + x + x2 + (1 | trial)
Df     AIC     BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
model2   4 19.4031 21.3428 -5.7016  11.4031
model2b  5  4.5587  6.9833  2.7206  -5.4413 16.844      1  4.057e-05 ***


improvement of confidence interval for estimate of slope, by using random effects:

> confint(model1)
2.5 %    97.5 %
(Intercept)  1.859806  3.275332
x           -3.459469 -1.116473
> confint(model2)
2.5 %    97.5 %
.sig01       0.2224486  1.493417
.sigma       0.1893400  0.487630
(Intercept)  1.7294434  3.405695
x           -2.8145888 -1.761353


So, this provides some numbers and lots of juggling with models...

...yet a simple visual peek on your data already provides this insight, as Peter Flom said

• Note that both a power law and exponential function provide good fits as well for this data:$$\text{model4 <- nls(y ~ a[trial] + b*exp(c*x), start=list(a=c(0,0,0), b=3, c=-1))}$$ $$\text{model5 <- nls(y ~ a[trial] + b*x^c, start=list(a=c(0,0,0), b=3, c=0.5))}$$With this little amount of data, and short range of variation in the x parameter (such that power law and exponential function are approximately linear), you can fit any func. It would be important to select a model, based on theoretic arguments (which is based on experience, previous observation), rather than this single experiment. – Sextus Empiricus Oct 19 '17 at 8:23