Comparision of multiple curves in R using linear mixed models

Question

first of all, I am sorry if a similar question already exits, in that case, please direct me to it!

I've got a repeated measurement of a value over 55-time points (720 min) in triplicate for two conditions (faint lines in the plot below). Now, I want to assess whether the mean curves (mean across all three replicates for each timepoint, bold lines in plot below) between the conditions show significantly different trajectories in R.

Example of the data structure is here:

     Timepoint  Time Value     N condition rep   concat_label
       <dbl> <dbl> <dbl> <dbl> <chr>     <chr> <chr>       
 1         0     0  269.   157 A         R1    A_R1        
 2         0     0  261.   123 B         R1    B_R1        
 3         1    10  426.   162 A         R1    A_R1        
 4         1    10  290.   126 B         R1    B_R1        
 5         2    20  632.   172 A         R1    A_R1        
 6         2    20  320.   126 B         R1    B_R1        
 7         3    30  753.   178 A         R1    A_R1        
 8         3    30  339.   130 B         R1    B_R1        
 9         4    40  790.   178 A         R1    A_R1        
10         4    40  350.   132 B         R1    B_R1  

In order to do so, I applied growth curve analysis and followed this tutorial (http://www.danmirman.org/gca). With that, I tried to derive a linear mixed model using the lme4 package and the lmer function. In my eyes, the model should have the structure:

lmer(Value ~ Time*Condition + (Time | rep)) to take the variability of the replicates into account. In order to get the right repression of the curves, I am fitting the linear model onto a fourth-order polynom using poly(Time, 4)) and can then convert the t-values into p-vals (I need to report a p-val although I know that they are debatable).

Here is an example code to perform the linear mixed model and generate a list of p-values from the t-values.

    model <- lmer(Value ~ poly(Time, 4)*condition + 
                  (poly(Time, 4) | rep),  
                control = lmerControl(optimizer="bobyqa"),
                data=data, REML=F)

coef(summary(model))

coefs2 <- data.frame(coef(summary(model))) 
coefs2$p <- format.pval(2*(1-pnorm(abs(coefs2$t.value))), digits=2)
coefs2

This returns values like:

> coefs2
                                Estimate Std..Error    t.value       p
(Intercept)                     792.4871  35.297235  22.451819 < 2e-16
poly(Timepoint, 4)1            1499.0161  97.490831  15.375971 < 2e-16
poly(Timepoint, 4)2            -664.3474  84.327948  -7.878140 3.3e-15
poly(Timepoint, 4)3             740.4694 112.671814   6.571914 5.0e-11
poly(Timepoint, 4)4            -586.8091  86.774484  -6.762462 1.4e-11
conditionB                     -357.6736   6.163614 -58.029855 < 2e-16
poly(Timepoint, 4)1:conditionB -667.3512 111.967603  -5.960217 2.5e-09
poly(Timepoint, 4)2:conditionB  542.1113 111.967603   4.841680 1.3e-06
poly(Timepoint, 4)3:conditionB -491.2822 111.967603  -4.387718 1.1e-05
poly(Timepoint, 4)4:conditionB  371.6393 111.967603   3.319168   9e-04

which I believe indicate that there is a significant difference in the trajectory of the example lines as indicated in "conditionB" (p < 2e-16). However, I am not sure whether this is the right model. Therefore, could someone give me some insights into whether this is the right model or how I might have to tweak in order to answer my specific question?

In addition to that, in a more advanced scenario, I would ideally like to compare multiple curves (e.g. conditions A vs B vs C vs D) to see which ones are significantly different (like in ANOVA situation but for non-linear data). Does someone have an idea on how to compare multiple curves using linear mixed models or even other types of regression analysis?

Answer 1

I have found this Cross Validated answer useful to "compare nonlinear fits". Here is how I understood the idea.

Let's say you have three groups of observations which are identified by the levels g1,g2, g3 of a factor grp. One model, let's call it M3, would be built by fitting the same nonlinear model to each of the three groups (i.e. you get "three fits in one"), the other model M0 would be a single fit to the pooled data where you ignore the group membership. Then you check with the anova() function whether M0 and M3 are "significantly different". If yes, then at least one of the curves is "different" from the rest. If not, then "grouping has no effect", M0 is just as good as M3, so any difference between the separately fitted curves "could be due to chance alone", they are not "significantly different".

This procedure uses a somewhat underdocumented feature of R's nls() function that lets you fit "indexed formulae". Basically, if your model is $$y(x) = \mathrm{MyFunc}(x; a, b)$$ where $a,b$ are the parameters, then the "indexed" formula would look like

y ~ MyFunc(x, a[grp], b[grp])

i.e. the parameters are "indexed" by the factor grp representing the groups and the result will be as if nls() had performed three fits on the three data subsets for which grp==g1, grp==g2, grp==g3, respectively.

I have used this trick on an enzyme inhibition dataset that contained Michaelis-Menten kinetics data points for an enzyme. The three groups were labeled "no inhibitor", and "low" and "high" concentrations of a competitive inhibitor. The question was whether the inhibitor has an effect.

I am not sure whether there's an easy way in R to do this with e.g. orthogonal polynomials, but maybe this idea will inspire you to try it out :-)