Pryseley Assam
2006-Mar-07 13:45 UTC
[R] lme and gls : accessing values from correlation structure and variance functions
Dear R-users
I am relatively new to R, i hope my many novice questions are welcome.
I have problems accessing some objects (specifically the random effects,
correlation structure and variance function) from an object of class gls and
lme.
I used the following models:
yah <- gls (outcome~ -1 + as.factor(Trial):as.factor(endpoint)+
as.factor(Trial):as.factor(endpoint):trt, data=datt[datt$Trial<4,],
correlation = corSymm(form=~1|as.factor(Trial)/as.factor(subject)),
weights=varIdent(form=~1|endpoint))
bm <- lme (outcome~ -1 + as.factor(endpoint)+ as.factor(endpoint):trt,
data=datt[datt$Trial<4,],
random=~-1 + as.factor(endpoint) + as.factor(endpoint):trt|as.factor(Trial),
correlation = corSymm(form=~1|as.factor(Trial)/as.factor(subject)),
weights=varIdent(form=~1|endpoint))
When i print the object "bm" i get the following output:
----------------------------------------------------------------------------------------------------------
> bm
Linear mixed-effects model fit by REML
Data: datt[datt$Trial < 4, ]
Log-restricted-likelihood: -52.23147
Fixed: outcome ~ -1 + as.factor(endpoint) + as.factor(endpoint):trt
as.factor(endpoint)-1 as.factor(endpoint)1 as.factor(endpoint)-1:trt
-3.663087 -1.772427 -3.661823
as.factor(endpoint)1:trt
-3.209671
Random effects:
Formula: ~-1 + as.factor(endpoint) + as.factor(endpoint):trt | as.factor(Trial)
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
as.factor(endpoint)-1 2.05744327 as.()-1 as.()1 a.()-1:
as.factor(endpoint)1 0.08400874 -0.976
as.factor(endpoint)-1:trt 1.90318009 0.975 -0.967
as.factor(endpoint)1:trt 3.25432832 -0.992 0.982 -0.972
Residual 6.48819860
Correlation Structure: General
Formula: ~1 | as.factor(Trial)/as.factor(subject)
Parameter estimate(s):
Correlation:
1
2 0.812
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | endpoint
Parameter estimates:
1 -1
1.000000 1.764878
Number of Observations: 18
Number of Groups: 3
----------------------------------------------------------------------------------------------------
Can somebody tell me how to access the values in blue above?
Also, when i tried accessing these values i obtained the following
bm$modelStruct
corStruct parameters:
[1] -1.394879
varStruct parameters:
[1] 0.5680815
What do these values represent.
Thanks in advance..
Pryseley
sample data:
RowNames Trial subject VISUAL0 TRT VISUAL24 VISUAL52 TREAT outcome endpoint
trt
4 1 1003 65 4 65 55 2 0 1 1
8 1 1007 67 1 64 68 2 -3 1 -1
12 2 1110 59 4 53 42 2 -6 1 1
14 2 1111 64 1 72 65 2 8 1 -1
16 2 1112 39 1 37 37 2 -2 1 -1
18 2 1115 59 4 54 58 2 -5 1 1
24 3 1806 46 4 27 24 2 -19 1 1
26 3 1813 31 4 33 48 2 2 1 1
28 3 1815 64 1 67 64 2 3 1 -1
4 1 1003 65 4 65 55 2 -10 -1 1
8 1 1007 67 1 64 68 2 1 -1 -1
12 2 1110 59 4 53 42 2 -17 -1 1
14 2 1111 64 1 72 65 2 1 -1 -1
16 2 1112 39 1 37 37 2 -2 -1 -1
18 2 1115 59 4 54 58 2 -1 -1 1
24 3 1806 46 4 27 24 2 -22 -1 1
26 3 1813 31 4 33 48 2 17 -1 1
28 3 1815 64 1 67 64 2 0 -1 -1
---------------------------------
[[alternative HTML version deleted]]
Spencer Graves
2006-Mar-11 18:19 UTC
[R] lme and gls : accessing values from correlation structure and variance functions
You ask, "how to access the values in blue". I can't see what you
carefully highlighted in blue, because this listserve often removes
formatting.
Have you consulted Pinheiro and Bates (2000) Mixed-Effects Models in
S and S-Plus (Springer)? I have learned a lot from that book, both
about how to use the nlme package and about mixed-effects modeling more
generally, and I would expect that you could find there answers to this
question and many others you will want answered in the future.
Beyond that, I checked the documentation for "gls" and
"lme" and
found they produced objects of class "gls" and "lme",
respectively. I
then entered 'methods(class="gls")' and
'methods(class="lme")' and got
the following:
> methods(class="lme")
[1] ACF.lme* anova.lme augPred.lme*
[4] BIC.lme* coef.lme* comparePred.lme*
[7] fitted.lme* fixef.lme* formula.lme*
[10] getData.lme* getGroups.lme* getGroupsFormula.lme*
[13] getResponse.lme* getVarCov.lme* intervals.lme*
[16] logLik.lme* pairs.lme* plot.lme
[19] predict.lme* print.anova.lme* print.lme*
[22] qqnorm.lme* ranef.lme* residuals.lme*
[25] simulate.lme summary.lme* update.lme*
[28] VarCorr.lme* Variogram.lme* vcov.lme*
Non-visible functions are asterisked
> methods(class="gls")
[1] ACF.gls* anova.gls* augPred.gls*
[4] BIC.gls* coef.gls* comparePred.gls*
[7] fitted.gls* formula.gls* getData.gls*
[10] getGroups.gls* getGroupsFormula.gls* getResponse.gls*
[13] getVarCov.gls* intervals.gls* logLik.gls*
[16] plot.gls* predict.gls* print.gls*
[19] qqnorm.gls* residuals.gls* summary.gls*
[22] update.gls* Variogram.gls* vcov.gls*
Non-visible functions are asterisked
>
This suggests to me that you might like to read '?VarCorr.lme'.
hope this helps.
spencer graves
Pryseley Assam wrote:
> Dear R-users
>
> I am relatively new to R, i hope my many novice questions are welcome.
>
> I have problems accessing some objects (specifically the random effects,
correlation structure and variance function) from an object of class gls and
lme.
>
> I used the following models:
>
>
> yah <- gls (outcome~ -1 + as.factor(Trial):as.factor(endpoint)+
as.factor(Trial):as.factor(endpoint):trt, data=datt[datt$Trial<4,],
> correlation = corSymm(form=~1|as.factor(Trial)/as.factor(subject)),
weights=varIdent(form=~1|endpoint))
>
>
> bm <- lme (outcome~ -1 + as.factor(endpoint)+ as.factor(endpoint):trt,
data=datt[datt$Trial<4,],
> random=~-1 + as.factor(endpoint) +
as.factor(endpoint):trt|as.factor(Trial),
> correlation = corSymm(form=~1|as.factor(Trial)/as.factor(subject)),
weights=varIdent(form=~1|endpoint))
>
> When i print the object "bm" i get the following output:
>
----------------------------------------------------------------------------------------------------------
> > bm
> Linear mixed-effects model fit by REML
> Data: datt[datt$Trial < 4, ]
> Log-restricted-likelihood: -52.23147
> Fixed: outcome ~ -1 + as.factor(endpoint) + as.factor(endpoint):trt
> as.factor(endpoint)-1 as.factor(endpoint)1
as.factor(endpoint)-1:trt
> -3.663087 -1.772427
-3.661823
> as.factor(endpoint)1:trt
> -3.209671
>
> Random effects:
> Formula: ~-1 + as.factor(endpoint) + as.factor(endpoint):trt |
as.factor(Trial)
> Structure: General positive-definite, Log-Cholesky parametrization
> StdDev Corr
> as.factor(endpoint)-1 2.05744327 as.()-1 as.()1 a.()-1:
> as.factor(endpoint)1 0.08400874 -0.976
> as.factor(endpoint)-1:trt 1.90318009 0.975 -0.967
> as.factor(endpoint)1:trt 3.25432832 -0.992 0.982 -0.972
> Residual 6.48819860
>
> Correlation Structure: General
> Formula: ~1 | as.factor(Trial)/as.factor(subject)
> Parameter estimate(s):
> Correlation:
> 1
> 2 0.812
>
> Variance function:
> Structure: Different standard deviations per stratum
> Formula: ~1 | endpoint
> Parameter estimates:
> 1 -1
> 1.000000 1.764878
> Number of Observations: 18
> Number of Groups: 3
>
----------------------------------------------------------------------------------------------------
>
> Can somebody tell me how to access the values in blue above?
>
> Also, when i tried accessing these values i obtained the following
>
> bm$modelStruct
> corStruct parameters:
> [1] -1.394879
> varStruct parameters:
> [1] 0.5680815
>
> What do these values represent.
>
> Thanks in advance..
> Pryseley
>
>
> sample data:
>
> RowNames Trial subject VISUAL0 TRT VISUAL24 VISUAL52 TREAT outcome
endpoint trt
> 4 1 1003 65 4 65 55 2 0 1 1
> 8 1 1007 67 1 64 68 2 -3 1 -1
> 12 2 1110 59 4 53 42 2 -6 1 1
> 14 2 1111 64 1 72 65 2 8 1 -1
> 16 2 1112 39 1 37 37 2 -2 1 -1
> 18 2 1115 59 4 54 58 2 -5 1 1
> 24 3 1806 46 4 27 24 2 -19 1 1
> 26 3 1813 31 4 33 48 2 2 1 1
> 28 3 1815 64 1 67 64 2 3 1 -1
> 4 1 1003 65 4 65 55 2 -10 -1 1
> 8 1 1007 67 1 64 68 2 1 -1 -1
> 12 2 1110 59 4 53 42 2 -17 -1 1
> 14 2 1111 64 1 72 65 2 1 -1 -1
> 16 2 1112 39 1 37 37 2 -2 -1 -1
> 18 2 1115 59 4 54 58 2 -1 -1 1
> 24 3 1806 46 4 27 24 2 -22 -1 1
> 26 3 1813 31 4 33 48 2 17 -1 1
> 28 3 1815 64 1 67 64 2 0 -1 -1
>
>
>
> ---------------------------------
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> 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