1
$\begingroup$

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)

enter image description here

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.

$\endgroup$
6
  • 1
    $\begingroup$ summary(lm(net_growth ~ dilution, yourData)) $\endgroup$
    – pogibas
    Commented Oct 18, 2017 at 9:48
  • $\begingroup$ 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. $\endgroup$ Commented Oct 18, 2017 at 9:52
  • $\begingroup$ "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. $\endgroup$
    – Roland
    Commented Oct 18, 2017 at 11:48
  • $\begingroup$ @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. $\endgroup$ Commented Oct 19, 2017 at 8:07
  • 1
    $\begingroup$ @MartijnWeterings I know. My comment is partly for the benefit of others who find this question through its title. $\endgroup$
    – Roland
    Commented Oct 19, 2017 at 8:11

3 Answers 3

4
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ And I would wonder why there seems to be a hint of the plotted values being bounded below at zero. $\endgroup$
    – mdewey
    Commented Oct 18, 2017 at 12:37
  • 1
    $\begingroup$ 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. $\endgroup$
    – Roland
    Commented Oct 19, 2017 at 6:34
  • $\begingroup$ 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? $\endgroup$
    – whuber
    Commented Feb 6, 2020 at 16:13
3
$\begingroup$

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))

plot with abline

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)

ggplot2 geom_smooth lm

$\endgroup$
5
  • 1
    $\begingroup$ Thanks a lot. So there is no need to check any assumption at this step ? $\endgroup$
    – Justine Courboules
    Commented Oct 18, 2017 at 10:32
  • 2
    $\begingroup$ @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. $\endgroup$ Commented Oct 18, 2017 at 10:37
  • $\begingroup$ 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. $\endgroup$
    – Eric Fail
    Commented Oct 18, 2017 at 10:47
  • $\begingroup$ @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? $\endgroup$
    – Justine Courboules
    Commented Oct 18, 2017 at 11:03
  • $\begingroup$ @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. $\endgroup$ Commented Oct 18, 2017 at 12:53
2
$\begingroup$

While others have already provided correct answers, I wish to provide some unsolicited advice.

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.

demonstration of different fits

# 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[function which(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)
legend(x=0.8, y=2.5, 
       legend=c("linear lm","linear lmer", "quadratic lmer"), 
       col=c(2,3,4), lwd=1)

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

comparison of adding quadratic term

> 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

$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Oct 19, 2017 at 8:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.