Your example is entirely too complicated for me to parse in the time
available, but I have a few questions that I hope might help:
First, have you examined str(fit.lme) plus all the other help pages
listed under "See Also" in the "lme" help page, especially
"lmeObject"?
With luck, this may answer your questions.
Second, are all your random effects nested, or are some crossed? If
the model is all standard nesting with no complicated inhomogeneity of
variance or correlation structure, you may be working too hard: The
'nlme' package is designed to make standard nested mixed-model analysis
as easily as Doug Bates could make it with a general tool in the S
language paradigm. Functions like 'pdBlocked' are relatively
inefficient attempts to adapt the 'nlme' paradigm to support crossed
random effects, etc. If you have both crossed and nested random
effects, you might try Bates' current project, 'lmer' associated
with
the 'lme4' package. Unfortunately, that package currently has fewer
helper functions for things like confidence intervals. (It's Bates'
leading edge development effort.)
Third, can generate a much simpler, self-contained problem that
exhibits the same features as your more complicated example? If your
email had contained a few lines of R code that someone like me could
copy and paste into R and get what you got, it would increase your
chances of getting a more informative reply quicker. Without it, I'm
left to guessing.
Fourth, if you want to know, "which is more correct", have you
considered Monte Carlo? The nlme package includes a function
simulate.lme, which might help you. The lme4 package includes mcmcsamp.
Hope this helps.
Spencer Graves
Kellie J. Archer, Ph.D. wrote:> I am trying to use lme to fit a mixed effects model to get the same
> results as when using the following SAS code:
>
> proc mixed;
> class refseqid probeid probeno end;
> model expression=end logpgc / ddfm=satterth;
> random probeno probeid / subject=refseqid type=cs;
> lsmeans end / diff cl; run;
>
> There are 3 genes (refseqid) which is the large grouping factor, with
> 2 probeids nested within each refseqid, and 16 probenos nested within
> each of the probeids.
>
> I have specified in the SAS Proc Mixed procedure that the
> variance-covariance structure is to be compound symmetric. Therefore,
> the variance-covariance matrix is a block diagonal matrix of the form,
>
> V_1 0 0
> 0 V_2 0
> 0 0 V3
>
> where each V_i represents a RefSeqID. Moreover, for each V_i the
> structure within the block is
>
> v_{11} v{12}
> v_{21} v{22}
>
> where v_{11} and v_{22} are different probeids nested within the
> refseqid, and so are correlated. The structure of
> both v_{11} and v_{22} are compound symmetric, and v_{12} and v{21}
> contain a constant for all elements of the matrix.
>
> I have tried to reproduce this using lme, but it is unclear from the
> documentation (and Pinheiro & Bates text) how the pdBlocked and
> compound symmetric structure can be combined.
>
> fit.lme<-lme(expression~End+logpgc,random=list(RefSeqID=pdBlocked(list
>
(~1,~ProbeID-1),pdClass="pdSymm")),data=dataset,correlation=corCompSym
> m(form=~1|RefSeqID/ProbeID/ProbeNo))
>
>
> The point estimates are essentially the same comparing R and SAS for
> the fixed effects, but the 95% confidence intervals are much shorter
> using lme(). In order to find the difference in the algorithms used by
> SAS and R I tried to extract the variance-covariance matrix to look at
> its structure. I used the getVarCov() command, but it tells me that
> this function is not available for nested structures. Is there another
> way to extract the variance-covariance structure for nested models?
> Does anyone know how I could get the var-cov structure above using
> lme?
>
>
> Kellie J. Archer, Ph.D.
> Assistant Professor, Department of Biostatistics
> Fellow, Center for the Study of Biological Complexity
> Virginia Commonwealth University
> 1101 East Marshall St., B1-066
> Richmond, VA 23298-0032
> phone: (804) 827-2039
> fax: (804) 828-8900
> e-mail: kjarcher at vcu.edu
> website: www.people.vcu.edu/~kjarcher
>
> ______________________________________________
> R-help at 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