Re: [R] lme vs. SAS proc mixed. Point estimates and SEs are the same, DFs are different
John Sorkin wrote:
R 2.3
Windows XP
I am trying to understand lme. My aim is to run a random effects regression
in which the intercept and jweek are random effects. I am comparing output
from SAS PROC MIXED with output from R. The point estimates and the SEs are
the same, however the DFs and the p values are different. I am clearly doing
something wrong in my R code. I would appreciate any suggestions of how I can
change the R code to get the same DFs as are provided by SAS.
This has been hashed over a number of times before. In short:
1) You're not necessarily doing anything wrong
2) SAS PROC MIXED is not necessarily doing it right
3) lme() is _definitely_ not doing it right in some cases
4) both work reasonably in large sample cases (but beware that this is
not equivalent to having many observation points)
SAS has an implementation of the method by Kenward and Rogers, which
could be the most reliable general DF-calculation method around (I don't
trust their Satterthwaite option, though). Getting this or equivalent
into lme() has been on the wish list for a while, but it is not a
trivial thing to do.
SAS code:
proc mixed data=lipids2;
model ldl=jweek/solution;
random int jweek/type=un subject=patient;
where lastvisit ge 4;
run;
SAS output:
Solution for Fixed Effects
Standard
Effect Estimate Error DFt ValuePr |t|
Intercept 113.48 7.4539 25 15.22 .0001
jweek -1.7164 0.5153 24 -3.33 0.0028
Type 3 Tests of Fixed Effects
Num Den
Effect DF DFF ValuePr F
jweek 1 24 11.090.0028
R code:
LesNew3 - groupedData(LDL~jweek | Patient, data=as.data.frame(LesData3),
FUN=mean)
fit3- lme(LDL~jweek, data=LesNew3[LesNew3[,lastvisit]=4,],
random=~1+jweek)
summary(fit3)
R output:
Random effects:
Formula: ~1 + jweek | Patient
Structure: General positive-definite, Log-Cholesky parametrization
Fixed effects: LDL ~ jweek
Value Std.Error DF t-value p-value
(Intercept) 113.47957 7.453921 65 15.224144 0.
jweek-1.71643 0.515361 65 -3.330535 0.0014
John Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
Confidentiality Statement:
This email message, including any attachments, is for the so...{{dropped}}
__
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
and provide commented, minimal, self-contained, reproducible code.
__
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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] lme vs. SAS proc mixed. Point estimates and SEs are the same, DFs are different
Peter Dalgaard wrote:
John Sorkin wrote:
R 2.3
Windows XP
I am trying to understand lme. My aim is to run a random effects regression
in which the intercept and jweek are random effects. I am comparing output
from SAS PROC MIXED with output from R. The point estimates and the SEs are
the same, however the DFs and the p values are different. I am clearly doing
something wrong in my R code. I would appreciate any suggestions of how I
can change the R code to get the same DFs as are provided by SAS.
This has been hashed over a number of times before. In short:
1) You're not necessarily doing anything wrong
2) SAS PROC MIXED is not necessarily doing it right
3) lme() is _definitely_ not doing it right in some cases
4) both work reasonably in large sample cases (but beware that this is
not equivalent to having many observation points)
SAS has an implementation of the method by Kenward and Rogers, which
could be the most reliable general DF-calculation method around (I don't
trust their Satterthwaite option, though). Getting this or equivalent
into lme() has been on the wish list for a while, but it is not a
trivial thing to do.
Forgot to say: All DF-based corrections are wrong if you have
non-normally distributed data (they depend on the 3rd and 4th moment of
the error distribution(s)), although they can be useful as warning signs
even in those cases. I also forgot to point to the simulate.lme()
function which can simulate the LR statistics directly.
SAS code:
proc mixed data=lipids2;
model ldl=jweek/solution;
random int jweek/type=un subject=patient;
where lastvisit ge 4;
run;
SAS output:
Solution for Fixed Effects
Standard
Effect Estimate Error DFt ValuePr |t|
Intercept 113.48 7.4539 25 15.22 .0001
jweek -1.7164 0.5153 24 -3.33 0.0028
Type 3 Tests of Fixed Effects
Num Den
Effect DF DFF ValuePr F
jweek 1 24 11.090.0028
R code:
LesNew3 - groupedData(LDL~jweek | Patient, data=as.data.frame(LesData3),
FUN=mean)
fit3- lme(LDL~jweek, data=LesNew3[LesNew3[,lastvisit]=4,],
random=~1+jweek)
summary(fit3)
R output:
Random effects:
Formula: ~1 + jweek | Patient
Structure: General positive-definite, Log-Cholesky parametrization
Fixed effects: LDL ~ jweek
Value Std.Error DF t-value p-value
(Intercept) 113.47957 7.453921 65 15.224144 0.
jweek-1.71643 0.515361 65 -3.330535 0.0014
John Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
Confidentiality Statement:
This email message, including any attachments, is for the so...{{dropped}}
__
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
and provide commented, minimal, self-contained, reproducible code.
__
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
and provide commented, minimal, self-contained, reproducible code.
__
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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] lme vs. SAS proc mixed. Point estimates and SEs are the same, DFs are different
In addition to Peter's comments, the following link summarizes the issue
as well. This is a direct response to the SAS/lmer DF issue.
https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html
Harold
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Peter Dalgaard
Sent: Tuesday, June 05, 2007 4:44 AM
To: John Sorkin
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] lme vs. SAS proc mixed. Point estimates and
SEs are the same, DFs are different
John Sorkin wrote:
R 2.3
Windows XP
I am trying to understand lme. My aim is to run a random
effects regression in which the intercept and jweek are
random effects. I am comparing output from SAS PROC MIXED
with output from R. The point estimates and the SEs are the
same, however the DFs and the p values are different. I am
clearly doing something wrong in my R code. I would
appreciate any suggestions of how I can change the R code to
get the same DFs as are provided by SAS.
This has been hashed over a number of times before. In short:
1) You're not necessarily doing anything wrong
2) SAS PROC MIXED is not necessarily doing it right
3) lme() is _definitely_ not doing it right in some cases
4) both work reasonably in large sample cases (but beware
that this is not equivalent to having many observation points)
SAS has an implementation of the method by Kenward and
Rogers, which could be the most reliable general
DF-calculation method around (I don't trust their
Satterthwaite option, though). Getting this or equivalent
into lme() has been on the wish list for a while, but it is
not a trivial thing to do.
SAS code:
proc mixed data=lipids2;
model ldl=jweek/solution;
random int jweek/type=un subject=patient;
where lastvisit ge 4;
run;
SAS output:
Solution for Fixed Effects
Standard
Effect Estimate Error DFt ValuePr |t|
Intercept 113.48 7.4539 25 15.22 .0001
jweek -1.7164 0.5153 24 -3.33 0.0028
Type 3 Tests of Fixed Effects
Num Den
Effect DF DFF ValuePr F
jweek 1 24 11.090.0028
R code:
LesNew3 - groupedData(LDL~jweek | Patient,
data=as.data.frame(LesData3), FUN=mean)
fit3- lme(LDL~jweek,
data=LesNew3[LesNew3[,lastvisit]=4,], random=~1+jweek)
summary(fit3)
R output:
Random effects:
Formula: ~1 + jweek | Patient
Structure: General positive-definite, Log-Cholesky parametrization
Fixed effects: LDL ~ jweek
Value Std.Error DF t-value p-value
(Intercept) 113.47957 7.453921 65 15.224144 0.
jweek-1.71643 0.515361 65 -3.330535 0.0014
John Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
Confidentiality Statement:
This email message, including any attachments, is for the
so...{{dropped}}
__
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
and provide commented, minimal, self-contained, reproducible code.
__
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
and provide commented, minimal, self-contained, reproducible code.
__
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
and provide commented, minimal, self-contained, reproducible code.
[R] lme vs. SAS proc mixed. Point estimates and SEs are the same, DFs are different
R 2.3
Windows XP
I am trying to understand lme. My aim is to run a random effects regression in
which the intercept and jweek are random effects. I am comparing output from
SAS PROC MIXED with output from R. The point estimates and the SEs are the
same, however the DFs and the p values are different. I am clearly doing
something wrong in my R code. I would appreciate any suggestions of how I can
change the R code to get the same DFs as are provided by SAS.
SAS code:
proc mixed data=lipids2;
model ldl=jweek/solution;
random int jweek/type=un subject=patient;
where lastvisit ge 4;
run;
SAS output:
Solution for Fixed Effects
Standard
Effect Estimate Error DFt ValuePr |t|
Intercept 113.48 7.4539 25 15.22 .0001
jweek -1.7164 0.5153 24 -3.33 0.0028
Type 3 Tests of Fixed Effects
Num Den
Effect DF DFF ValuePr F
jweek 1 24 11.090.0028
R code:
LesNew3 - groupedData(LDL~jweek | Patient, data=as.data.frame(LesData3),
FUN=mean)
fit3- lme(LDL~jweek, data=LesNew3[LesNew3[,lastvisit]=4,],
random=~1+jweek)
summary(fit3)
R output:
Random effects:
Formula: ~1 + jweek | Patient
Structure: General positive-definite, Log-Cholesky parametrization
Fixed effects: LDL ~ jweek
Value Std.Error DF t-value p-value
(Intercept) 113.47957 7.453921 65 15.224144 0.
jweek-1.71643 0.515361 65 -3.330535 0.0014
John Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
Confidentiality Statement:
This email message, including any attachments, is for the so...{{dropped}}
__
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
and provide commented, minimal, self-contained, reproducible code.
Re: [R] lme and SAS Proc mixed
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
00 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: [EMAIL PROTECTED]
website: www.people.vcu.edu/~kjarcher
__
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
__
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
[R] lme and SAS Proc mixed
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
00 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: [EMAIL PROTECTED]
website: www.people.vcu.edu/~kjarcher
__
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
