# Linear regression with repeated measures in R

I was unable to figure out how to perform linear regression in R in for a repeated measure design. In a previous question (still unanswered) it was suggested to me to not use lm but rather to use mixed models. I used lm in the following way:

lm.velocity_vs_Velocity_response <- lm(Velocity_response~Velocity*Subject, data=mydata)


(more details on the dataset can be found at the link above)

However I was not able to find on the internet any example with R code showing how to perform a linear regression analysis.

What I want is on one hand a plot of the data with the line fitting the data, and on the other hand the $R^2$ value along with the p-value for the test of significance for the model.

Is there anyone who can provide some suggestions? Any R code example could be of great help.

Edit
According to the suggestion I received so far, the solution to my analyze my data in order to understand if there is a linear relation between the two variables Velocity_response (deriving from the questionnaire) and Velocity (deriving from the performance) should be this:

library(nlme)
summary(lme(Velocity_response ~ Velocity*Subject, data=scrd, random= ~1|Subject))


The result of summary gives this:

    > summary(lme(Velocity_response ~ Velocity*Subject, data=scrd, random= ~1|Subject))
Linear mixed-effects model fit by REML
Data: scrd
AIC      BIC   logLik
104.2542 126.1603 -30.1271

Random effects:
Formula: ~1 | Subject
(Intercept) Residual
StdDev:    2.833804 2.125353

Fixed effects: Velocity_response ~ Velocity * Subject
Value Std.Error DF    t-value p-value
(Intercept)               -26.99558  25.82249 20 -1.0454288  0.3083
Velocity                   24.52675  19.28159 20  1.2720292  0.2180
SubjectSubject10           21.69377  27.18904  0  0.7978865     NaN
SubjectSubject11           11.31468  33.51749  0  0.3375754     NaN
SubjectSubject13           52.45966  53.96342  0  0.9721337     NaN
SubjectSubject2           -14.90571  34.16940  0 -0.4362299     NaN
SubjectSubject3            26.65853  29.41574  0  0.9062674     NaN
SubjectSubject6            37.28252  50.06033  0  0.7447517     NaN
SubjectSubject7            12.66581  26.58159  0  0.4764880     NaN
SubjectSubject8            14.28029  31.88142  0  0.4479188     NaN
SubjectSubject9             5.65504  34.54357  0  0.1637076     NaN
Velocity:SubjectSubject10 -11.89464  21.07070 20 -0.5645111  0.5787
Velocity:SubjectSubject11  -5.22544  27.68192 20 -0.1887672  0.8522
Velocity:SubjectSubject13 -41.06777  44.43318 20 -0.9242591  0.3664
Velocity:SubjectSubject2   11.53397  25.41780 20  0.4537754  0.6549
Velocity:SubjectSubject3  -19.47392  23.26966 20 -0.8368804  0.4125
Velocity:SubjectSubject6  -29.60138  41.47500 20 -0.7137162  0.4836
Velocity:SubjectSubject7   -6.85539  19.92271 20 -0.3440992  0.7344
Velocity:SubjectSubject8  -12.51390  22.54724 20 -0.5550080  0.5850
Velocity:SubjectSubject9   -2.22888  27.49938 20 -0.0810519  0.9362
Correlation:
(Intr) Velcty SbjS10 SbjS11 SbjS13 SbjcS2 SbjcS3 SbjcS6 SbjcS7 SbjcS8 SbjcS9 V:SS10 V:SS11 V:SS13 Vl:SS2 Vl:SS3
Velocity                  -0.993
SubjectSubject10          -0.950  0.943
SubjectSubject11          -0.770  0.765  0.732
SubjectSubject13          -0.479  0.475  0.454  0.369
SubjectSubject2           -0.756  0.751  0.718  0.582  0.362
SubjectSubject3           -0.878  0.872  0.834  0.676  0.420  0.663
SubjectSubject6           -0.516  0.512  0.490  0.397  0.247  0.390  0.453
SubjectSubject7           -0.971  0.965  0.923  0.748  0.465  0.734  0.853  0.501
SubjectSubject8           -0.810  0.804  0.769  0.624  0.388  0.612  0.711  0.418  0.787
SubjectSubject9           -0.748  0.742  0.710  0.576  0.358  0.565  0.656  0.386  0.726  0.605
Velocity:SubjectSubject10  0.909 -0.915 -0.981 -0.700 -0.435 -0.687 -0.798 -0.469 -0.883 -0.736 -0.679
Velocity:SubjectSubject11  0.692 -0.697 -0.657 -0.986 -0.331 -0.523 -0.607 -0.357 -0.672 -0.560 -0.517  0.637
Velocity:SubjectSubject13  0.431 -0.434 -0.409 -0.332 -0.996 -0.326 -0.378 -0.222 -0.419 -0.349 -0.322  0.397  0.302
Velocity:SubjectSubject2   0.753 -0.759 -0.715 -0.580 -0.360 -0.992 -0.661 -0.389 -0.732 -0.610 -0.563  0.694  0.528  0.329
Velocity:SubjectSubject3   0.823 -0.829 -0.782 -0.634 -0.394 -0.622 -0.984 -0.424 -0.799 -0.667 -0.615  0.758  0.577  0.360  0.629
Velocity:SubjectSubject6   0.462 -0.465 -0.438 -0.356 -0.221 -0.349 -0.405 -0.995 -0.449 -0.374 -0.345  0.425  0.324  0.202  0.353  0.385
Velocity:SubjectSubject7   0.961 -0.968 -0.913 -0.740 -0.460 -0.726 -0.844 -0.496 -0.986 -0.778 -0.718  0.886  0.674  0.420  0.734  0.802
Velocity:SubjectSubject8   0.849 -0.855 -0.807 -0.654 -0.406 -0.642 -0.746 -0.438 -0.825 -0.988 -0.635  0.783  0.596  0.371  0.649  0.709
Velocity:SubjectSubject9   0.696 -0.701 -0.661 -0.536 -0.333 -0.526 -0.611 -0.359 -0.676 -0.564 -0.990  0.642  0.488  0.304  0.532  0.581
Vl:SS6 Vl:SS7 Vl:SS8
Velocity
SubjectSubject10
SubjectSubject11
SubjectSubject13
SubjectSubject2
SubjectSubject3
SubjectSubject6
SubjectSubject7
SubjectSubject8
SubjectSubject9
Velocity:SubjectSubject10
Velocity:SubjectSubject11
Velocity:SubjectSubject13
Velocity:SubjectSubject2
Velocity:SubjectSubject3
Velocity:SubjectSubject6
Velocity:SubjectSubject7   0.450
Velocity:SubjectSubject8   0.398  0.828
Velocity:SubjectSubject9   0.326  0.679  0.600

Standardized Within-Group Residuals:
Min          Q1         Med          Q3         Max
-1.47194581 -0.46509026 -0.05537193  0.39069634  1.89436646

Number of Observations: 40
Number of Groups: 10
Warning message:
In pt(q, df, lower.tail, log.p) : NaNs produced
>


Now, I do not understand where I can get the R^2 and the corresponding p-values indicating me wether there is a linear relationship between the two variables or not, nor I have understood how my data can be plotted with the line fitting the regression.

Can anyone be so kind to enlighten me? I really need your help guys...

-
"Mixed Effects Models and Extensions in Ecology with R" by Zuur et al. is a nice introduction to linear mixed effects models, which focusses less on theory and more on application of the methodology. –  Roland Sep 3 '12 at 9:26
Dear Roland, I believe that that book is useful, but I am rather searching for something on line...do you have any webpage to suggest? –  Luca Sep 3 '12 at 10:43
As I said in your earlier post, lm() has a plot associated with it. So, if your model is M1 you can use plot(M1). –  Peter Flom Sep 3 '12 at 10:51
dear @PeterFlom yes, but you also told me to avoid using lm for repeated measures design. So, my question is if I have to use lm for analyzing my data or another function. Any suggestion? –  Luca Sep 3 '12 at 11:02
Like I said, look into multi-level models. In R, you can look at the nlmepackage. Also, search this site for the topic, there is a lot written about it here. –  Peter Flom Sep 3 '12 at 11:27
show 1 more comment

What you do really depends on the goals of the analysis. I'm not certain exactly what the goals of your analysis are, but I'll go through several examples, and hopefully one of them will be applicable to your situation.

Case 1: One quantitive variable measured twice

Let's say that you ran a human subject study in which you had participants take a stats test twice and you wanted to find out if the average scores on the second measurement were different from the first measurement (to determine whether learning occurred). If scores test1 and test2 are stored in data frame d, You could do this entirely using the lm() function, as in:

mod <- lm(test2 - test1 ~ 1, data = d)
summary(mod)


The test of the intercept is the test of the difference between test1 and test2. Note that you will have no delta-R^2 for the difference between test1 and test2 -- instead, your measure of effect size should be something like cohen's d.

Case 2a: One quantitative variable measured twice, one dichotomous variable, measured totally between subjects

Let's say that we have the same study design, but we want to know if different rates of learning occurred for men and women. So, we have one quantitative variable (test performance) that is measured twice, and one dichotomous variable, measured once. Assuming test1, test2, and gender are all contained in data frame d, We could also test this model only using lm(), as in:

mod <- lm(test2 - test1 ~ gender, data = d)
summary(mod)
lm.sumSquares(mod) # lm.sumSquares() is located in the lmSupport package, and gives the change in R^2 due to the between-subjects part of the model


Assuming gender is centered (i.e., coded, for example, male = -.5 and female = +.5), the intercept in this model is the test of the difference between test 1 and test 2, averaged across males and females. The coefficient for gender is the interaction between time and gender. To get the effect of gender, averaged across time, you would have to do:

mod <- lm(rowMeans(cbind(test2, test1)) ~ gender, data = d)
summary(mod)


Case 2b: One quantitative variable measured twice, one quantitative variable, only measured once

Let's assume that again we have one quantitative variable measured twice and one quantitative variable measured once. So, for example, let's say we had a measure of baseline interest in statistics and we wanted to determine whether people who had higher levels of baseline interest learned more from time 1 to time 2. We'd first have to center interest, as in:

d$interestc <- d$interest - mean(d$interest)  Assuming that test1, test2, and interestc are all in data frame d, this question could then be tested very similarly to Case 1a: mod <- lm(test2 - test1 ~ interestc, data = d) summary(mod) lm.sumSquares(mod)  Once again, the intercept in this model tests whether, averaged across interest, test scores changed from time 1 to time 2. However, this interpretation only holds when interest is centered. The coefficient for interest is whether the effect of time depends on baseline interest. We could get the effect of interest, averaged across time, by averaging together test1 and test 2, as above, and testing the effect of interest on this composite variable. Case 2c: One quantitative variable measured twice, one categorical variable, only measured once Let's assume that your between-subjects variable was a category, measured only once. So, for example, let's assume that you were interested in whether people of different races (White vs Asian vs Black vs Hispanic) had different amounts of learning from time 1 to time 2. Assuming test1, test2, and race are in data frame d, you would first need to contrast code race. This could be done using planned orthogonal contrasts, dummy codes, or using effects codes, depending on specific hypotheses / questions you want to test (I recommend looking at lm.setContrasts() if you're looking for a helper function to do this). Assuming the race variable is already contrast coded, you would use lm() very similarly to the above two cases, as in: mod <- lm(test2 - test1 ~ race, data = d) summary(mod) lm.sumSquares(mod)  Assuming the race contrasts are centered, the intercept in this model is, once again, the "main effect" of time. The coefficients for the race contrasts are the interactions between those contrasts and time. To obtain the omnibus effects of race, use the following code: Anova(mod, type = 3)  Case 3: One quantitative variable measured 3 times (i.e., a three-level within-subjects manipulation) Let's assume that you added a third point of measurement to the design from case one. So, your participants took a stats test three times instead of twice. Here you have a couple choices, depending on whether you want an omnibus test of the differences between time points (sometimes you don't). For example, let's say that your main hypothesis is that test scores will linearly increase from time 1 to time 3. Assuming that test1, test2, and test3 are in data frame d, this hypothesis could be tested by first creating the following composite: d$lin <- d[, paste("test", sep = "", 1:3)] %*% c(-1, 0, 1)


Then you would test whether an intercept only model using lin as the dependent variable has an intercept that is different from 0, as in:

mod <- lm(lin ~ 1, data = d)
summary(mod)


This will give you your test of whether stats scores were increasing over time. You can, of course, create other types of custom difference scores, depending on your particular hypotheses.

If you care about omnibus tests of significance, you need to use the Anova() function from the car package. The specific implementation is a little convoluted. Basically, you specify which variables are within-subjects and which are between-subjects using lm(). You then create the within-subjects portion of the model (i.e., specify which of test1, test2, and test3 were measured first, second, and third) and then pass that model to Anova() by creating a data frame called idata. Using my hypothetical example:

mod <- lm(cbind(test1, test2, test3) ~ 1, data = d) # No between-subjects portion of the model
idata <- data.frame(time = c("time1", "time2", "time3")) # Specify the within-subjects portion of the model
mod.A <- Anova(mod, idata = idata, idesign = ~time, type = 3) # Gives multivariate tests.  For univariate tests, add multivariate = FALSE
summary(mod.A)


The idesign statement tells Anova to include the time variable (composed of test1, test2, and test3) in the model. This code will give you your omnibus tests of the effects of time on test scores.

Case 4: One quantitative variable measured 3 times, one between-subjects quantitative variable

This case is a simple extension of Case 3. As above, if you merely care about 1 degree of freedom tests, you can simply create a custom difference score with your within-subjects variable. So, assuming that test1, test2, test3, and interest are all in data frame d, and assuming that we are interested in the linear effects of time on test scores (and how those effects of time vary with baseline interest), you would do the following:

d\$lin <- d[, paste("test", sep = "", 1:3)] %*% c(-1, 0, 1)


Then, do the following (with mean-centered interest):

mod <- lm(lin ~ interestc, data = d)
summary(mod)
lm.sumSquares(mod)


If you want omnibus tests, do the following:

mod <- lm(cbind(test1, test2, test3) ~ interest, data = d) # We now have a between-subjects portion of the model
idata <- data.frame(time = c("time1", "time2", "time3"))
mod.A <- Anova(mod, idata = idata, idesign = ~time * interest, type = 3) # The idesign statement assumes that we're interested in the interaction between time and interest
summary(mod.A)


Other cases: I will omit these for brevity, but they are simple extensions of what I've already described.

Please note that the (univariate) omnibus tests of time where time has more than 2 levels all assume sphericity. This assumption becomes pretty untenable as you increase the number of levels. If you have quite a few points of measurement in your design (say, 4+) I strongly recommend you use something like multilevel modeling and move to a package that is specialized for this technique (such as nlme or lme4.

Hope this helps!

-
Dear Patrick @user1188407 , thanks a lot you have been very kind in providing such answer. Unfortunately my case probably fits in what you wrote in the last sentences...so I would need a R code example to understand how to treat my data. Indeed if you look at the design of my experiment described in a previous post stackoverflow.com/questions/12182373/… you can see that I have a variable measured 4 time (i.e. the velocity measured in 4 conditions) –  Luca Sep 3 '12 at 22:33
and I want to find if there is a linear relation with a variable (velocity_response) expressing the perceived velocity in the four conditions. So each participant underwent 4 conditions and then evaluated the perception of those conditions. I want to know if the performance is related with the perception... –  Luca Sep 3 '12 at 22:35
Well if you want to use a multilevel modeling solution, you can draw on a lot of different online resources. To start with, you should take a look at the nlme package and this vignette. The vignette is slightly out of date (2002), I found it useful when I was learning about multi-level modeling. Finally, you can check out the book published by the makers of the nlme package. –  Patrick S. Forscher Sep 4 '12 at 2:47
Dear Patrick @user1188407 thanks. I studied the multilevel models, and I arrived to this formula to analyze my data: lme(Velocity_response ~ Velocity*Subject, data=scrd, random= ~1|Subject) Can you please confirm me that this formula is correct for the analysis I want to perform on my data? However, I do not understand how can I get R^2 and p-values, nor how to plot the graphic with the points and the line fitting the regression. Could you please help me? I am not a statician... –  Luca Sep 4 '12 at 11:23
The formula seems correct to me based on my understanding of your study (and assuming you've formatted your data in person period format). You would get your p-values by saving the results of your analysis into an object (as I do in my examples) and getting a summary of that object. However, because of the differences between multilevel models and traditional regression (e.g., in effect size metrics -- there's no standard metric in multilevel models) I strongly suggest you read more about this technique before using it. It seems like the other users have recommended several good options. –  Patrick S. Forscher Sep 4 '12 at 13:21