Linear regression with repeated measures in R

Question

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

Answer 1

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!