5
$\begingroup$

I have a data-set with two independent (continuous) variables: temperature (1, 3, 6 degrees) and nutrient concentration (~ 9-1400) and one dependent variable (length), with replicates (three independent replicates at each condition). I'd like to analyze the interactive effects of nutrients and temperature on size.

When I plot the data, it is clear (visually) that length increases with nutrient concentration (to a point), and temperature doesn't play a big rule in influencing length.

enter image description here

However, I would like to use some statistical test to back up my observations.

require(lme4) 
require(nlme)
require(lmerTest)
require(gam)

#this is the data I have
temperature <- c(1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6)

nutrient <- c(41.922282,41.922282,41.922282,37.23794,37.23794,37.23794,31.662541,31.662541,31.662541,279.720746,279.720746,279.720746,248.465109,248.465109,248.465109,211.264016,211.264016,211.264016,27.946784,27.946784,27.946784,24.824046,24.824046,24.824046,21.1073,21.1073,21.1073,55.899183,55.899183,55.899183,49.65308,49.65308,49.65308,42.218842,42.218842,42.218842,97.838316,97.838316,97.838316,86.905988,86.905988,86.905988,73.89411,73.89411,73.89411,10.479381,10.479381,10.479381,9.308428,9.308428,9.308428,7.914737,7.914737,7.914737,1404.251692,1404.251692,1404.251692,1247.342413,1247.342413,1247.342413,1060.585803,1060.585803,1060.585803,139.790093,139.790093,139.790093,124.170127,124.170127,124.170127,105.578928,105.578928,105.578928)

replicate <- c("1C_41.922282","1C_41.922282","1C_41.922282","3C_37.23794","3C_37.23794","3C_37.23794","6C_31.662541","6C_31.662541","6C_31.662541","1C_279.720746","1C_279.720746","1C_279.720746","3C_248.465109","3C_248.465109","3C_248.465109","6C_211.264016","6C_211.264016","6C_211.264016","1C_27.946784","1C_27.946784","1C_27.946784","3C_24.824046","3C_24.824046","3C_24.824046","6C_21.1073","6C_21.1073","6C_21.1073","1C_55.899183","1C_55.899183","1C_55.899183","3C_49.65308","3C_49.65308","3C_49.65308","6C_42.218842","6C_42.218842","6C_42.218842","1C_97.838316","1C_97.838316","1C_97.838316","3C_86.905988","3C_86.905988","3C_86.905988","6C_73.89411","6C_73.89411","6C_73.89411","1C_10.479381","1C_10.479381","1C_10.479381","3C_9.308428","3C_9.308428","3C_9.308428","6C_7.914737","6C_7.914737","6C_7.914737","1C_1404.251692","1C_1404.251692","1C_1404.251692","3C_1247.342413","3C_1247.342413","3C_1247.342413","6C_1060.585803","6C_1060.585803","6C_1060.585803","1C_139.790093","1C_139.790093","1C_139.790093","3C_124.170127","3C_124.170127","3C_124.170127","6C_105.578928","6C_105.578928","6C_105.578928")

length <- c(0.284222,0.271812,0.287842,0.266703,0.325212,0.323167,0.368914,0.307848,0.331279,0.349361,0.344158,0.379752,0.418207,0.398789,0.397851,0.481935,0.46838,0.447341,0.291471,0.38784,0.355353,0.353436,0.40762,0.321866,0.284687,0.26343,0.281308,0.361157,0.367518,0.328645,0.390822,0.372086,0.366396,0.357013,0.388808,0.440506,0.351289,0.348172,0.345575,0.35433,0.363403,0.332073,0.34037,0.315966,0.351829,0.207838,0.227385,0.183385,0.198436,0.217075,0.270751,0.28564,0.228815,0.212524,0.410496,0.415918,0.416817,0.406967,0.38017,0.417732,0.453175,0.502706,0.477136,0.371708,0.344421,0.366723,0.398991,0.393513,0.442445,0.414689,0.442346,0.446943)

mydata <-data.frame(nutrient, temperature, replicate, length)

#to keep things simple, I am using a linear mixed effects model, and using my "replicates" as random effects. 
#I also try different transformations to handle the non-linearity in my data

linear <- lmer (length ~ nutrient*temperature +(1|replicate),  data= mydata)
log <- lmer (length ~ log(nutrient)*temperature +(1|replicate),  data= mydata)
sqrt <- lmer (length ~ sqrt(nutrient)*temperature +(1|replicate),  data= mydata)

#then I do an AIC test to see which of these linear transformations is best 

AIC(linear, log, sqrt)

#this shows that the log transformation is the best model. 

Then I take the log model, and run it with nutrient alone, temperature alone, nutrient*temperature, and nutrient+temperature effects: 

log_interactive <- lmer (length ~ log(nutrient)*temperature +(1|replicate),  data= mydata)
log_additive <- lmer(length ~ log(nutrient)+temperature +(1|replicate),  data= mydata)
log_temperature <- lmer (length ~ temperature +(1|replicate),  data= mydata)
log_nutrient <- lmer (length ~ log(nutrient) +(1|replicate),  data= mydata)

AIC(log_interactive, log_additive, log_temperature, log_nutrient)
#this shows that the best model is the nutrient alone one. 

My interpretation of this is that nutrients play a role in influencing size, temperature does not influence size, and there is no temperature-nutrient interaction in influencing size.

I have a very weak statistical background and I am not sure that this is the correct approach to take to analyze this data and generate p-values for the different effects. Any help would be appreciated. Would it be better to do non-linear statistics on this ?

Thanks in advance!

$\endgroup$
  • 4
    $\begingroup$ Your data will be easier to see if you use a logarithmic scale for the nutrient axis, and it may be more usefully analysed. See this question and answers for more information: stats.stackexchange.com/questions/27951/… $\endgroup$ – Michael Lew Jan 25 at 20:22
2
$\begingroup$

Often with data that look like this, the goal is to determine the x value beyond which there is no further (statistical) increase in the y variable. And also to determine the plateau y value. This might be done with linear-plateau or quadratic-plateau models. But eyeballing your data, a Cate-Nelson approach might be more useful to determine these values. This approach basically separates the low-x-low-y values from the high-x-high-y values.

As a first approach, I might look at the data for each temperature separately. I might estimate the critical x value and plateau value for each of these separately, and then determine if the temperatures have a statistical effect by comparing the 95% confidence intervals of these statistics. It looks to me that Temperature 6 has a higher plateau than Temperature 1.

For the purposes of this answer, I'm just going to examine the data for Temperature 6.


First, let's try a linear plateau model on the Temperature 6 data.

temperature <- c(1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6,1,1,1,3,3,3,6,6,6)

nutrient <- c(41.922282,41.922282,41.922282,37.23794,37.23794,37.23794,31.662541,31.662541,31.662541,279.720746,279.720746,279.720746,248.465109,248.465109,248.465109,211.264016,211.264016,211.264016,27.946784,27.946784,27.946784,24.824046,24.824046,24.824046,21.1073,21.1073,21.1073,55.899183,55.899183,55.899183,49.65308,49.65308,49.65308,42.218842,42.218842,42.218842,97.838316,97.838316,97.838316,86.905988,86.905988,86.905988,73.89411,73.89411,73.89411,10.479381,10.479381,10.479381,9.308428,9.308428,9.308428,7.914737,7.914737,7.914737,1404.251692,1404.251692,1404.251692,1247.342413,1247.342413,1247.342413,1060.585803,1060.585803,1060.585803,139.790093,139.790093,139.790093,124.170127,124.170127,124.170127,105.578928,105.578928,105.578928)

replicate <- c("1C_41.922282","1C_41.922282","1C_41.922282","3C_37.23794","3C_37.23794","3C_37.23794","6C_31.662541","6C_31.662541","6C_31.662541","1C_279.720746","1C_279.720746","1C_279.720746","3C_248.465109","3C_248.465109","3C_248.465109","6C_211.264016","6C_211.264016","6C_211.264016","1C_27.946784","1C_27.946784","1C_27.946784","3C_24.824046","3C_24.824046","3C_24.824046","6C_21.1073","6C_21.1073","6C_21.1073","1C_55.899183","1C_55.899183","1C_55.899183","3C_49.65308","3C_49.65308","3C_49.65308","6C_42.218842","6C_42.218842","6C_42.218842","1C_97.838316","1C_97.838316","1C_97.838316","3C_86.905988","3C_86.905988","3C_86.905988","6C_73.89411","6C_73.89411","6C_73.89411","1C_10.479381","1C_10.479381","1C_10.479381","3C_9.308428","3C_9.308428","3C_9.308428","6C_7.914737","6C_7.914737","6C_7.914737","1C_1404.251692","1C_1404.251692","1C_1404.251692","3C_1247.342413","3C_1247.342413","3C_1247.342413","6C_1060.585803","6C_1060.585803","6C_1060.585803","1C_139.790093","1C_139.790093","1C_139.790093","3C_124.170127","3C_124.170127","3C_124.170127","6C_105.578928","6C_105.578928","6C_105.578928")

length <- c(0.284222,0.271812,0.287842,0.266703,0.325212,0.323167,0.368914,0.307848,0.331279,0.349361,0.344158,0.379752,0.418207,0.398789,0.397851,0.481935,0.46838,0.447341,0.291471,0.38784,0.355353,0.353436,0.40762,0.321866,0.284687,0.26343,0.281308,0.361157,0.367518,0.328645,0.390822,0.372086,0.366396,0.357013,0.388808,0.440506,0.351289,0.348172,0.345575,0.35433,0.363403,0.332073,0.34037,0.315966,0.351829,0.207838,0.227385,0.183385,0.198436,0.217075,0.270751,0.28564,0.228815,0.212524,0.410496,0.415918,0.416817,0.406967,0.38017,0.417732,0.453175,0.502706,0.477136,0.371708,0.344421,0.366723,0.398991,0.393513,0.442445,0.414689,0.442346,0.446943)

mydata <-data.frame(nutrient, temperature, replicate, length)

T6 = mydata[mydata$temperature==6,]

plot(length ~ nutrient, data=T6, pch=16, col="firebrick")

Plot

###  Guess some reasonable initial values for parameters

a.ini     = 0
b.ini     = 1000
clx.ini   = 200

###  Define linear plateau function

linplat = function(x, a, b, clx)
          {ifelse(x < clx, a + b * x,
                           a + b * clx)}

###  Find best fit parameters

model = nls(length ~ linplat(nutrient, a, b, clx),
            data = T6,
            start = list(a   = a.ini,
                         b   = b.ini,
                         clx = clx.ini),
             trace = FALSE,
             nls.control(maxiter = 1000))

summary(model)

   ### Parameters:
   ###      Estimate Std. Error t value Pr(>|t|)    
   ### a   2.607e-01  1.711e-02  15.231 8.01e-13 ***
   ### b   1.618e-03  2.972e-04   5.446 2.11e-05 ***
   ### clx 1.305e+02  1.957e+01   6.665 1.35e-06 ***

plateau = 2.607e-01 + 1.618e-03 * 1.305e+02
plateau

   ### 0.472

So the critical x value is about 130 and the plateau is 0.472. This suggests that above a nutrient value of 130, there is no further rise in length with increasing nutrient value (for the Temperature 6 data).

With the caveat that I am the author of the page, additional code to determine confidence intervals for the parameters, p value for the overall model, and pseudo r-square for the model can be found at rcompanion.org/handbook/I_11.html.

This model is easy enough to plot, but I'll use a convenience function from the rcompanion package.

require(rcompanion)

plotPredy(data  = T6,
          x     = nutrient,
          y     = length,
          model = model,
          xlab  = "Nutrient",
          ylab  = "Length for T 6")

Plot

This model might be okay for these data. But in my mind, the data below nutrient = 200 don't really support a linear rise in length in relation to nutrient. The residuals also reveal that the model doesn't quite fit the trend of the data in this region.

plot(predict(model), residuals(model))

hist(residuals(model), col="darkgray")

We could try a quadratic-plateau model. As it turns out, it's not much different.

Plot

Another idea is to use a Cate-Nelson approach. The goal here is to find a critical x value and a critical y that separates the low-x-low-y values from the high-x-high-y values. There are statistical approaches that can be employed. But if there are cutoff values of x or y that may be meaningful, this analysis can be accomplished by simply eyeballing these values. (An example might be, we are interested in getting at least 85% of maximum yield of our crop, so we want to determine the nutrient level in the soil (x) that gives us (usually) at least 85% of maximum yield (y).)

I won't get into the specifics of analysis, but there are some resources here: rcompanion.org/rcompanion/h_02.html, (with the caveat that I am the author of that page).

The following plots one possible solution.

Plot

[Note the points in this plot are slightly jittered from the original data.]

It turns out that this solution perfectly separates the data. That is, all data points fall into Quadrants II and IV.

The Critical x value is 158, and the Critical y is 0.447. That is a nutrient value greater than 158 is likely to yield a length of greater than 0.447. And vice-versa.

sum((T6$nutrient < 158) & (T6$length < 0.447142))

   ### 18

sum((T6$nutrient > 158) & (T6$length > 0.447142))

   ### 6

sum((T6$nutrient < 158) & (T6$length > 0.447142))

   ### 0

sum((T6$nutrient > 158) & (T6$length < 0.447142))

   ### 0
$\endgroup$
  • $\begingroup$ Thank you for your reply! So, just to make sure I'm understanding this correctly. The output here shows that length increases to a maximum of ~0.447. Nutrients increase cause an increase in length, but not after nutrient levels are above 158. Would you advise that I run this separately for the 3 different temperatures, and compare the outputs ? $\endgroup$ – Loay Jabre Jan 26 at 17:22
  • $\begingroup$ Yes, your description of the results describe the logic of the linear-plateau and quadratic-plateau models, though the numbers don't match. Note that the Cate-Nelson approach doesn't really ensure there's any "increase"; it just says "there are a bunch of low-x-low-y values and a bunch of high-x-high-y values, and this is one way to separate them". (Cate-Nelson can also be conducted in a more "statistical" manner based on the sums of squares of the groups, but I didn't get in to that here.) $\endgroup$ – Sal Mangiafico Jan 27 at 17:43
  • $\begingroup$ I think I would look at the different temperatures separately. If you can get confidence intervals for the critical x and plateau levels, and if all these overlap for the different temperatures, you could say they're not different, and pool all the data. That would be my approach, but I'm not committed to this approach. $\endgroup$ – Sal Mangiafico Jan 27 at 17:46
  • $\begingroup$ BTW, I think fitting a logarithmic curve as @rvl suggests, also makes sense. One consideration is what you think happens in reality biologically: Should length plateau at high nutrient levels? Or should length continue to increase a small amount at high levels? A 2nd consideration is if you are interested in determining some kind of critical x level above which there is no or little increase in y... If you fit a smooth curve, you might look into a Mitscherlich-Bray model, which will look similar to the simple logarithmic model, but has a theoretical maximum y value. $\endgroup$ – Sal Mangiafico Jan 28 at 11:23
2
$\begingroup$

If you do as was suggested in a comment by Michael Lew and use log(nutrient) in the plot, you see a reasonably linear trend for each temperature:

lattice::xyplot(length ~ log(nutrient) | temperature, data = mydata)

scatterplot with log(nutrient)

Accordingly, I tried the following model, which fits quite well:

mylm = lm(length ~ factor(temperature) * log(nutrient), data = mydata)

The fitted trends can be visualized as follows:

library(emmeans)
emmip(mylm, temperature ~ nutrient, at = list(nutrient = c(10,20,40,80,160,320,640,1280)))

The fitted lines This appears as linear trends because the values I chose go up exponentially. If we look at the same plot with a linear scale...

 pdat = .Last.value$data
lattice::xyplot(yvar ~ nutrient, groups = ~temperature, data = pdat, type="l")

Plot of fitted curves

You can compare the slopes of the lines in the preceding plot as follows:

> emtrends(mylm, pairwise ~ temperature, var = "log(nutrient)")
$emtrends
 temperature log(nutrient).trend      SE df lower.CL upper.CL
           1              0.0350 0.00614 66   0.0228   0.0473
           3              0.0317 0.00614 66   0.0194   0.0439
           6              0.0523 0.00614 66   0.0400   0.0645

Confidence level used: 0.95 

$contrasts
 contrast estimate      SE df t.ratio p.value
 1 - 3     0.00333 0.00868 66  0.383  0.9223 
 1 - 6    -0.01728 0.00868 66 -1.992  0.1222 
 3 - 6    -0.02061 0.00868 66 -2.375  0.0527 

P value adjustment: tukey method for comparing a family of 3 estimates

I think this model fits the data better than broken lines or polynomials. It is parsimonious, and the interpretation is simple.

$\endgroup$
  • $\begingroup$ Another reason, BTW, for choosing a model with log(nutrient) is that the leverage of the data points is more equitably distributed. See the first plot, $\endgroup$ – Russ Lenth Jan 27 at 16:20
1
$\begingroup$

The first step usually taken for non-linear relationships like this is by adding a polynomial term. A linear relationship is y ~ x, but you can use a quadratic approximation by doing y ~ x + x ^ 2. This is like x interacting with itself: the relationship between x and y depends on the value of x. It is still a linear model, since you are adding up coefficients, but it is a very common "trick" to approximate nonlinear relationships while still using all of the benefits of the linear model.

Here's code how to do it with lme4 and then the plotted relationship:

library(lme4) 
library(lmerTest)

# linear nutrient and temperature main effects
# treating temperature as numeric
m1 <- lmer(length ~ nutrient + temperature + (1 | replicate), mydata)
summary(m1) # significant main effects of nutrient, but not temeprature

# looks like curvilinear relationship between nutrient and length
# simplest way to do this is quadratic: x + x ^ 2
# you can think of x ^ 2 as "x times x" or an interaction between x and itself
# poly() gives us polynomial, and we indicate 2 for squared
m2 <- lmer(length ~ poly(nutrient, 2) + temperature + (1 | replicate), mydata)
summary(m2) # significant quadratic relationship

# does this quadratic relationship depend on temperature?
m3 <- lmer(length ~ poly(nutrient, 2) * temperature + (1 | replicate), mydata)
summary(m3)
# yes, p = .033. we see the quadratic effect of nutrient and
#   linear of temperature depend on one another. what does this look like? 
#   let's get predicted values for the data and plot it

library(ggplot2)
library(emmeans)
pred <- emmeans(
  m3, # specify model 
  c("nutrient", "temperature"), # specify predictors
  # specify at what points you want to get predictions:
  at = list(
    nutrient = seq(8, 1000, length.out = 20),
    temperature = c(1, 3, 6)
  )
)

# make data frame for plotting
pred <- as.data.frame(pred)

# and plot
ggplot(pred, aes(x = nutrient, y = emmean, 
                 group = factor(temperature),color = factor(temperature))) +
  geom_line() +
  labs(y = "predicted length")

m3 has the quadratic relationship between nutrient and length that depends on temperature, and you can see that coefficient is significant. Just the fixed effect coefficient table from summary(m3):

Fixed effects:
                                Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                     0.333554   0.017941 17.999958  18.591 3.39e-13 ***
poly(nutrient, 2)1              0.472626   0.162941 17.999958   2.901  0.00953 ** 
poly(nutrient, 2)2             -0.077723   0.150528 17.999958  -0.516  0.61190    
temperature                     0.004344   0.004663 17.999958   0.932  0.36386    
poly(nutrient, 2)1:temperature -0.103052   0.053974 17.999958  -1.909  0.07230 .  
poly(nutrient, 2)2:temperature -0.109066   0.047116 17.999958  -2.315  0.03263 *  

And here's what the call to ggplot returns:

enter image description here

If you are interested in what the different predicted effects for replicate look like, you could also add (1 + replicate) to the call to emmeans and include the IDs that you want to look at.

$\endgroup$
  • 1
    $\begingroup$ You can also use ggeffects() to get the marginal effects and plot for the graph with a lot less code: library(ggeffects) to load the library and then ggpredict(m3, c("nutrient [all]", "temperature")) %>% plot() to get the plot. If you remove [all] the plot will look exactly like @MarkWhite's. $\endgroup$ – Erik Ruzek Jan 25 at 21:14
  • $\begingroup$ Thank you for the response! When I plot the raw data on top of this: ggplot(pred, aes(x = nutrient, y = emmean, group = factor(temperature),color = factor(temperature))) + geom_point()+ geom_point(data = mydata, aes(x = nutrient, y = length))+ geom_line() + labs(y = "predicted length") It looks like the models don't really fit the data very well. It looks like it's overestimating the difference in temperature. Is this right ? $\endgroup$ – Loay Jabre Jan 25 at 21:34
  • 2
    $\begingroup$ I don’t think a polynomial model will be very suitable for these data, because polynomials want very badly to go to infinity at each end of the scale, whereas the data plot looks like it levels off. A nonlinear model such as a one-compartment model (from pharmakinetics) seems like a better choice. $\endgroup$ – Russ Lenth Jan 25 at 22:42
  • $\begingroup$ @rvl sure, but in my experience comparing polynomial terms in linear models vs. something that theoretically fits the data better, the linear model is pretty robust to violations of assumptions in simulations. At least in terms of significance, which is what this question concerns. $\endgroup$ – Mark White Jan 25 at 23:04
  • 2
    $\begingroup$ I think replacing one unflexible functional form (a line) with another (a parabola) as "the first step" is bad science: let the data speak. $\endgroup$ – Alexis Jan 26 at 16:26

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.