Re: [R] LME & data with complicated random & correlational structures

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Wed 07 Dec 2005 - 13:32:52 EST

          Have you received any replies to this post? I haven't seen any, so I will attempt a few comments. First, I'm overwhelmed with the details in your discussion. I suggest that for each of your question, you try to think of an extremely simple example that would test your question, as suggested in the Posting Guide (www.R-project.org/posting-guide.html). If you have problems with that, please submit a question focused on that one thing.

          Have you looked at Pinheiro and Bates (2000) Mixed Effects Models in S and S-Plus (Springer)? If no, I suggest you take a hard look at this book. I was unable to get anything sensible from "lme" until I started working through this book. This is the primary reference for "lme", and for me at least, it was essential to understanding what I needed to do to get anything useful from "lme". You don't follow all the math in this book to get something useful out of it, because it includes many worked examples that should go a long way to answering many questions you might have about "lme".

  1. "if I set something in my model formula as a fixed effect, then it does not make sense to set it as a random effect as well? ...

 > Temperature ~ Stimulus-1, random=Subj|Subj/Epoch/Stimulus

          I perceive several problems here. Have you tried the following:

 > Temperature ~ Stimulus-1, random=~1|Subj/Epoch

 From your description, I'm guessing that Stimulus is a factor with levels like ("Baseline", "A", "B", "Recovery"). My way of thinking about these things is to try to write this as an algebraic model with parameters to estimate. "Temperature~Stimulus-1", ignoring the "random" argument for the moment, could be written as follows:

(1) Temperature = b["Baseline"]*I(Stimulus=="Baseline") + b["A"]*I(Stimulus=="A") + b["B"]*I(Stimulus=="B") + b["Recovery"]*I(Stimulus=="Recovery"),

          where the 4 "b" parameters are to be estimated by "iterative, reweighted generalized least squares" to maximize an appropriate likelihood function [and where I(...) = indicator function that is 1 if the (...) is TRUE and 0 otherwise].

          Now consider "random=~1|Subj"; ignore "epoch" for the moment. This adds one "random coefficient" to this model for each Sub. If you have 2 subjects, then the model becomes something like the following:

(2) Temperature = b["Baseline"]*I(Stimulus=="Baseline") + b["A"]*I(Stimulus=="A") + b["B"]*I(Stimulus=="B") + b["Recovery"]*I(Stimulus=="Recovery") + b.subj[1]*I(Subj==1) + b.subj[2]*I(Subj==2).

          However, we do NOT initially estimate the random coefficients b.subj[1] and b.subj[2]. Rather, we assume these coefficients "b.subj" are normally distributed with mean 0 and a variance "var.subj", and we want to estimate "var.subj". (We may later estimate the b.subj's

conditioned on the estimate of "var.subj", but that's a separate issue.) 
  We estimate "var.subj" using "iterative, reweighted generalized least 
squares" that roughly speaking "uncorrelates" or "whitens" the 
correlated residuals from model (1) and then minimizes the sum of squares of those "whitened" residuals. More detail is provided in Pinheiro and Bates.

          Now consider "random=~1|Subj/Epoch". This adds another random coefficient for each (Subject, Epoch) combination, which we also assume are normally distributed with mean 0 and variance "var.s.epoch". We then want to estimate the 4 fixed effect coefficients and the two variance coefficients simultaneously, using the same approach I just outlined.

        > 2- Is it possible to take a piecewise approach wrt the variance using lme(), such as modeling the variability of each subject first ... . When I try to set up the correlation structure, I run out of memory fast."

          You've hit on a great idea here: Consider the data for each subject separately. Plot it, think about the similarities and differences, and condense the data into, e.g., 1 or a few numbers per subject / epoch combination. Then use "lme" on this condensed data set. I'd start simple and add complexity later. I have on occasion tried to fit the most complicated model I could think of, only to find that I could have gotten most of the information from condensing it grossly and fitting much simpler models, and I wasted lots of time and effort trying to work with all the data at once. I suggest you start by aggregating it to the grossest level you think might answer your research questions. After you have sensible answers at that level, if you still have time and energy for this project, I then might try aggregate the data into more summaries with fewer observations in each summary.

          3. Is there a way to get rid of the serial dependency BEFORE running the model with LME(), such as initiating a corStruct before placing it in the model?

          Answer: Yes, and the primary way to do this would be to aggregate, like I just suggested. If you still want to do more with the correlation structure, study Pinheiro and Bates and try something with a moderately condensed version of what you have.

	  hope this helps.
	  spencer graves

Keith Chamberlain wrote:

> Dear List,
>
> This is my first post, and I'm a relatively new R user trying to work out a
> mixed effects model using lme() with random effects, and a correlation
> structure, and have looked over the archives, & R help on lme, corClasses, &
> etc extensively for clues. My programming experience is minimal (1 semester
> of C). My mentor, who has much more programming experience, but a comparable
> level of knowledge of R, has been stumped. What follows are 3 questions
> pertaining to an lme() model, one on the nested hierarcy, 1 on a strategy
> for a piecewise approach to the variance given I have ~24 hours of data
> (sampled at 32Hz, 1hr per subject), and one on the corStruct or how to get
> rid of serial dependencies before lme().
>
> I'm analyzing skin temperature continuously recorded at 32Hz in Baseline (10
> min), Testing (~5 min), and Recovery (20 min) epochs of a face recognition
> experiment. Stimuli are the same in Baseline and Recovery (portrait or
> landscape), and in testing, participants were tested on their recognition of
> a list b&w portraits presented just before testing started. On some of the
> portraits 'learned' the eyes were masked, and in others, the eyes were
> visible. In testing, the portraits have no masking but the stimuli in
> testing are labeled "Eyes" and "NoEyes". The data structure looks as
> follows:
>
> Subj/Epoch/Stimuli/Time/Temperature
> There are 8 subjects
>
> 9 epochs - 6 of which were just "instruction" blocks, and one "Learning"
> period. Wrt lme(), I figured out how to use subset too isolate just the
> Baseline, Learning, and Testing Epochs (and avoid epochs with only 1
> stimulus level, such as "instruction"). Data within each epoch are balanced
> wrt # trials, but not between epochs. Recovery has twice as many trials as
> Baseline, and Testing has about half. Time for each epoch is roughly that
> ratio too, although time in each trial differs.
>
> Stimuli are the same in Baseline & Recovery, but different in Testing,
> although there are 2 levels in each used epoch.
>
> Time & Temperature make up the time series, and Temperature is the dependent
> variable too stimulus.
>
> 1- are fixed effects and random effects discrete? That is, if I set
> something in my model formula as a fixed effect, then it does not make sense
> to set it as a random effect as well? The documentation (and posts) were not
> really clear on that point (not that the documentation technically 'should'
> be per say, just that I got confused).
>
> The nested hierarchy for what actually gets analyzed looks as follows:
> Subj/Epoch/Stimulus/Temperature
>
> Reasoning: there are several temperature samples recorded in each trial of
> Stimulus. Several stimuli in each Epoch, and all of the Epochs for one
> subject. Subject is random (theoretically) because of sampling in a
> population, Epoch would be fixed because all participants went through the
> same sequence of Epochs, but Stimulus varied randomly within an Epoch, which
> seems inconsistent when I apply it to the lme model as both a fixed and
> random effect.
>
> Temperature ~ Stimulus-1, random=Subj|Subj/Epoch/Stimulus
> Subset= Epoch=="Baseline" | Epoch=="Testing" | Epoch=="Recovery"
>
> I'm looking to correctly allocate error terms for between subjects (Subj)
> variability, and further delineate the within subject error between Epoch
> and Stimulus. The current model that I got to work (memory issues namely) is
> Temperature ~ Stimulus-1, random=Subj|Subj, which I decided to use to get
> the residuals to have the Subject variability accounted for and subtracted.
> Would a list of random structures work better? If so, is each item in the
> list structured just as the random formula? I haven't actually seen/found
> any examples of a list of random/nesting structures.
>
> 2- Is it possible to take a piecewise approach wrt the variance using lme(),
> such as modeling the variability of each subject first, then using
> further-nested terms in a model and the residuals from the previous? If so,
> what caveats exist for interpreting the variances?
>
> I'm not interpreting p-values at this point because of another issue. When I
> try to set up the correlation structure, I run out of memory fast. I've
> tried this on a mac G5, an HP Pavilion dv1000 (= Pentium 2.6GHz), and a
> Gateway with an AMD athalon 900MHz processors. Each system has 386M memory
> or more, one of which has 1G.
>
> 3- Is there a way to get rid of the serial dependency BEFORE running the
> model with LME(), such as initiating a corStruct before placing it in the
> model? I'm working with so much data that I'm fine with doing the process
> piecewise. An AR process was difficult because the residuals are not the
> same length as the data file that I started with. Serial dependencies still
> gota go, whether via the correlation term in lme() or some other method,
> because I'll soon be breaking up the variance into components via
> spectrum().
>
> So I might as well add a 4th. What's the term that gets me too data after
> AR() has done it's work? I'm thinking that resid() wasn't right but data
> that the data differ from their original length prior to an AR process may
> be how its done.
>
> Rgds,
> KeithC.
> Psych Undergrad, CU Boulder
> RE McNair Scholar
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves@pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Wed Dec 07 13:43:12 2005

This archive was generated by hypermail 2.1.8 : Wed 07 Dec 2005 - 20:23:43 EST