5
$\begingroup$

I have some doubts about the covariance structure in a multilevel model fitted in R (using the nlme package). I'm not an expert (just starting to learn statistics...), so I apologize if some of my questions seem evident. I've checked the previous posts and haven't found an answer.

I have data from an experiment in which we have registered physiological data from 30 subjects in 2 conditions (with 30 trials in each condition). These 30 trials are close in time, and we expected a higher correlation between closer trials, with the correlation decreasing as trials are further apart from each other. We are interested in the effect of condition, not in time effects. I think that the right way to analyze these data is to fit a multilevel model, in which TRIAL is a level 1 variable, SUBJECT is a level 2 variable, CONDITION is a fixed factor, and the DV is the physiological response (FR). The R command I'm using is:

lme(fixed= FR ~ CONDITION, data=mydata, random= ~ TIME | SUBJECT)

My (many) questions are both theoretical and practical:

  • Which covariance structure lme does use by default? Is it a problem not to use the most appropriate covariance structure?

  • I've read that an autoregressive covariance structure (AR1) refers to a constant variance at each time point and a weaker correlation as time points get further apart. My data only meets the second criterion. How can I know which covariance structure is right for my data? How important is it to the validity of the results?

  • I'm only interested in the CONDITION effect, which is significant when I'm not including TIME in the model. I'm not interesting in TIME effects, nor I want to use the model to make a prediction, but only to check the significance of the CONDITION effects. Is it correct if I drop out TIME and fit a model only with CONDITION?

Thank you for your help!

Improve this question
$\endgroup$
3
  • 4
    $\begingroup$ For question 1: in ?lme it gives: correlation an optional corStruct object describing the within-group correlation structure. ... Defaults to NULL, corresponding to no within-group correlations. So default is uncorrelated repeated measurements. $\endgroup$
    – Seth
    Jun 6, 2012 at 16:29
  • 1
    $\begingroup$ AR(1) with the correlation parameter rho less than 1 in absolute value is stationary. That means that the mean and variance do not change with time. The autocorrealtion function has an exponential decay declining from a positive value down to 0 as the lag increases when rho >0. If rho<0 it decays but alternates sign from negative for n=1 to positive for n=2 and back and forth negative for odd lags and positive for even ones. $\endgroup$
    – Michael R. Chernick
    Jun 6, 2012 at 19:39
  • $\begingroup$ There is a pretty good explanation of the use of different covariance structures within the lme() function at: rpsychologist.com/r-guide-longitudinal-lme-lmer $\endgroup$
    – hagenaue
    Jul 10, 2016 at 21:20

2 Answers 2

Reset to default
5
$\begingroup$

I'm not sure I can provide the kind of answer I'd like to, but I will try to throw out some pieces of information regarding your questions.

First, both @Seth and @gui11aume (+1 to each) have noted that lme() defaults to no within group correlations. The question is why, and whether that's likely to be a problem. I believe that the thinking is a properly specified multilevel model will account for the covariance amongst your observations such that the residuals are independent. That's why the function was coded to expect no correlations. That is, you may be OK.

Several of your questions concern the effect of having a misspecified variance/covariance structure (bearing in mind that this may not actually apply to you). The estimation of your betas should be unaffected by this, that is, they should be unbiased. However, the estimation of the variance of the sampling distributions will be inaccurate, that is, your p-values will be inaccurate. Moreover, I believe that you cannot say a-priori whether they will be too high or too low. If you are really concerned about these issues you can always use robust (a.k.a., 'sandwich') standard errors. These are typically thought about in the context of generalized linear models, but they can be used elsewhere. Check out the R package sandwich. Note that if they are not necessary, you could be at risk of increased type II errors.

The standard AR(1) variance/covariance structure does assume homoskedasticity, so far as I know. More restrictive, however, is that it assumes every observation was made at the appropriate time, and that all measurements are equally spaced in time. These assumptions usually don't hold, even in the most fortuitous situations, and as such, the AR(1) variance/covariance structure is dangerous to assume.

Remember that the proper specification of the model for the means is crucial. It is remotely possible that time is not relevant to the appropriate model of the mean, but it isn't very likely at all. Leaving TIME out of the model risks the omitted variable bias. Thus, dropping TIME is likely to yield both biased estimates of the means and invalid inferences. This is just not worth gambling on.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thank you @gui11aume and gung for your very useful answers (+1 to each). Concerning to the inclusion of TIME in the model, I understand and agree with your argumentations: it have to be in the model. I've tried to include it both as fixed and random factor, and the model is improved (-2LogLik comparison), although CONDITION is not a significant predictor. But I still keep doubts about the variance-covariance structure: I think my data are heterokedastic (unequal variances), so I should not specify AR1. But when I do it the model fits better. Is it a proof of anything? (continue) $\endgroup$
    – Mike_999
    Jun 8, 2012 at 7:50
  • $\begingroup$ On the other hand, there is an evident correlation between pairs of obervations, that decreases as time gets further apart. So I think my data are not completely unestructured. Is there any variance-covariance structure that accounts for this different variance at each time, but with correlations between the closer measures? Thank you! $\endgroup$
    – Mike_999
    Jun 8, 2012 at 7:55
  • 1
    $\begingroup$ I don't know why you think the data are within Subject heteroskedastic. It doesnt matter if the data are between subject heteroskedastic to fit the ar(1) model. I am not even sure how you detect within subject heteroskedasticity. In ar(1) is ok so long as the variance between all the measurements of one person are the same and you only have at most one pair of measurements between a subjects say 4th and 5th observation. Basically ar(1) is ok if your time intervals are equal length. $\endgroup$
    – Seth
    Jun 8, 2012 at 12:50
  • 1
    $\begingroup$ Mike, to notify someone (other than the question asker, you, or the answerer, me), you must put the 'at' symbol before their name. So, to notify Seth, you write @Seth. Next, I think what he's getting at is: are you finding that there are correlations amongst your observations, or among your residuals? Only the residuals matter. If you have a properly specified multilevel model, your residuals should be independent & homoskedastic, & if they are, then you are OK. $\endgroup$
    – gung - Reinstate Monica
    Jun 8, 2012 at 18:04
  • 1
    $\begingroup$ Mike, Say you have two time points and you plot something like a spaghetti plot where you have a line for each person and they basically start all together at time one and then spread out at time 2, I think that is what you are describing. This is not within subject heteroskedasticity, it is between subject heteroskedasticity. It often just means there is a high variation between subjects in outcome. But ar(1) can still work. Do like @gung suggests and check the residuals from a model without Time slopes, colored by ID and see if they display a linear (or quadratic) pattern. $\endgroup$
    – Seth
    Jun 8, 2012 at 18:15
5
$\begingroup$

You first point has been addressed by @Seth: the default is no within group correlation. Regarding your second point, I think it depends on the effect size of TIME. If the effect is minor compared to your other variables, like CONDITION and SUBJECT, the misfit should not be a problem. But if TIME is a major effect you will want your model to describe it well.

Now regarding your third question, if you drop TIME from the model by calling

lme(fixed= FR ~ CONDITION, data=mydata, random= ~ 1 | SUBJECT)

you remove the interaction term TIME*SUBJECT. So if TIME has the same effect for every subject, it's no big deal. However, if the response shows (or is expected to show) a particular behaviour for some people at some time (say some people get better and some get worse) then this will be "absorbed" in your other variables in ways that are difficult to predict (i.e. might make some terms significant and others non significant).

That you are not interested in TIME does not mean that you should leave it out of the model. Quite the contrary I would say. To give a heuristic example, in simple linear regression, you fit a slope and and intercept, but you are usually interested in the slope only. If you do not include the intercept term, you fit a a line through the origin. You can check by yourself that with this approach you will sometimes (for different clouds of points) reach the wrong conclusion: conclude that the slope is non null when it is and vice versa.

So if you have a doubt, keep TIME in but only focus one the estimates and p-values for CONDITION.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.