# How can I compute regression for several longitudinal data sets (thus, with auto-correlated error)?

My actual project is a bit complicated, but I'll explain by analogy (which I hope facilitates response):

I have 3 substances, say water, motor oil, and ethanol. For each substance, I have 5 samples in a beaker (total 15 beakers). I heat all the beakers on a hot-plate up to 70 degrees Celsius, and over the next hour, I measure the temperature of the fluid in each beaker at 5 minute intervals.

Newtonian cooling provides me with a good prediction about these temperature data, namely that the temperature of the fluid in each cup should follow an exponential distribution: y = a + e^(-kt) where a is room temperature.

I want to estimate the value k for each substance and test the hypothesis that k1 > k2 > k3 (1, 2, 3 corresponding to my three substances). The natural method of estimating k seems to be computing a non-linear regression on each substance's data, or possibly log-transforming all the data and then just computing a simple linear regression. However, there are problems.

Some questions:

1. Given the obvious autocorrelation in the longitudinal data (confirmed by my (P)ACF plots of course), must I compute an AR term and filter my data prior to computing the regression?
2. Assuming I compute this autoregression term, how do I compute it for five independent sets of data (the five beakers of a give substance)? I could average the five beakers together and then compute the regression, but this screws up my AR term (assuming I need one) and also throws off my estimate of the actual within-beaker variance from the model.
3. What completely wrong-headed assumption(s) have I worked in here...?
• Wouldn't you need a $b$ parameter in the function like $y = a + b \cdot e^{-kt}$? Where $b=70° - a$. – GaBorgulya Apr 10 '11 at 2:56
• Since your data has time trend it is natural that it will be autocorrelated. This does not mean however that you really need to correct for AR. Try fitting model without any special error variance assumption and then check for signs of autocorrelation. – mpiktas Apr 11 '11 at 8:35

As we have strong reasons to believe that the cooling will follow the $y(t) = a + e^{-kt}$ function for each beaker I would first check if this model fits the data well indeed.

If it does I wouldn't bother with analysing the autocorrelation at all, but focus on the estimation of $k_1$, $k_2$ and $k_3$, and testing the hypothesis about them.

To estimate $k_1$, $k_2$ and $k_3$ you need a non-linear model. Your idea of log transformation followed by linear modelling is best when the error (difference between the measured $y$ temperature and the one predicted by the formula) is proportional to the temperature. However, I suspect that the error will be primarily due to temperature measurement and thus normally distributed with the same variance for any temperature (you need to check this). If so, a non-linear model would be more appropriate.

A model using the above function will give you estimates for the parameters of the cooling of a single beaker, $a$ and $k$. We may however assume that $a$ should be the same for each beaker, $k$s should be similar for the same substance, and that the standard deviation ($\sigma$) is the same across all temperature measurements. These can be expressed in a model accounting for all the beakers in the same time (second index $j$ is beaker ID): $$y_j(t) = a + e^{-(k_i + \alpha_j)t} + \epsilon$$ where $\epsilon$ is normally distributed error of SD $\sigma$, $k_i$ is one of 3 mean $k$ values for substance $i$, $\alpha_j$ is normally distributed random deviation of a specific beaker from the $k_i$ substance mean, with a substance specific SD ($\sigma_{\alpha{}i}$). This is now a non-linear mixed effect model, that can be fitted using various software. After this you have the $k_i$ values and their standard errors.

The next question is how to test the hypothesis that $k_1 > k_2 > k_3$. It may be “cleaner” to formulate such a hypothesis in the Bayesian way. However you used the word test, so you probably want a significance test – but in order to do that you have to have a more specific alternative hypothesis (or family of hypotheses).

• Based on my fleeting understanding of dependent errors, not using an AR term is going to throw off the confidence intervals for my estimations, isn't it? Let's just say I want to build a bunch of 95% confidence intervals for each parameter and say "are they different?" (allowing for Bonferonni correction). Aren't the s.e.'s of k (as calculated by R or Matlab) way too conservative if my errors have dependency? – ben Apr 10 '11 at 2:46
• If the $y = a + e^{-(k_i + \alpha_j)t} + \epsilon$ model correctly describes the cooling the $\epsilon_{j,t}$ values should be independent. – GaBorgulya Apr 10 '11 at 2:58

If I understand your question correctly, you should be able to achieve what you want to do using a non-linear mixed-effects model. If you use R, you can use the nlme package. Basically as fixed factors you have a covariate (a) and a factor (substance or $i$ in $k_{i}$). You also have a random effect (individual measurements units or unitID). The good thing about nlme is that it also allows you to model the correlations in the residuals with e.g. an AR covariance structure.

edit: I always like to use a mixed-model when dealing with repeated measures. Still, if you don't want to include a random factor, you can model it with gnls in the same package. gnls still lets you select AR as covariance structure of he residuals.