[R-sig-ME] mcmcsamp
Andrew Robinson
A.Robinson at ms.unimelb.edu.au
Mon Dec 8 20:51:30 CET 2008
Hi Frank,
thanks ... but, I guess I should have asked you for a commentary as
well. Can you make our lives easier by identifying exactly what your
problem is?
Andrew
On Mon, Dec 08, 2008 at 10:24:05AM -0500, Frank Lawrence wrote:
> Hi Andrew:
>
> Sorry for not including the example at the outset.
>
> >sessionInfo()
> R version 2.8.0 (2008-10-20)
> i386-pc-mingw32
>
> locale:
> LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
> States.1252;LC_MONETARY=English_United
> States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
>
> attached base packages:
> [1] grid splines stats graphics grDevices datasets
> [7] tcltk utils methods base
>
> other attached packages:
> [1] MCMCpack_0.9-5 coda_0.13-3 statmod_1.3.8
> [4] polycor_0.7-6 sfsmisc_1.0-6 mvtnorm_0.9-2
> [7] xtable_1.5-4 prettyR_1.3-5 lme4_0.999375-27
> [10] Matrix_0.999375-16 effects_2.0-0 nnet_7.2-44
> [13] mvnormtest_0.1-6 xlsReadWrite_1.3.2 gmodels_2.14.1
> [16] gtools_2.5.0 latticeExtra_0.5-4 lattice_0.17-17
> [19] RColorBrewer_1.0-2 doBy_3.6 foreign_0.8-29
> [22] Design_2.1-2 survival_2.34-1 e1071_1.5-18
> [25] class_7.2-44 car_1.2-9 mitools_2.0
> [28] MASS_7.2-44 svSocket_0.9-5 TinnR_1.0.2
> [31] R2HTML_1.59 Hmisc_3.4-4
>
> ##artificial data
> > nn <- 1e2
>
> > mm <- seq(1,5,1)
>
> > cv <- matrix(data = rep(x = 0.3, times = 25), nc = 5, nr = 5)
>
> > diag(cv) <- 1
>
> > dat <- cbind.data.frame(id = seq(1,nn,1), mvrnorm(n = nn, m = mm, S = cv,
> emp = T))
>
> > names(dat)[2:6] <- paste('y',1:5,sep='')
>
> > d.t <- reshape(data = dat, varying = list(names(dat)[2:6]), v.names = 'y',
> times = seq(0,4,1), idvar = 'id', drop = NULL, dir = 'l')
>
> > m1 <- lmer(form = y ~ time + (1|id), data = d.t, fam = gaussian, R = F, na
> = na.exclude)
> > m1
> Linear mixed model fit by maximum likelihood
> Formula: y ~ time + (1 | id)
> Data: d.t
> AIC BIC logLik deviance REMLdev
> 1370 1387 -681 1362 1371
> Random effects:
> Groups Name Variance Std.Dev.
> id (Intercept) 0.409 0.640
> Residual 0.955 0.977
> Number of obs: 500, groups: id, 100
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 1.0000 0.0860 11.6
> time 1.0000 0.0309 32.4
>
> Correlation of Fixed Effects:
> (Intr)
> time -0.719
>
> > x <- mcmcsamp(obj = m1, n = 1e3)
>
> > str(x)
> Formal class 'merMCMC' [package "lme4"] with 9 slots
> ..@ Gp : int [1:2] 0 100
> ..@ ST : num [1, 1:1000] 0.655 0.529 0.448 0.401 0.406 ...
> ..@ call : language lmer(formula = y ~ time + (1 | id), data = d.t,
> REML = F, na.action = na.exclude)
> ..@ deviance: num [1:1000] 1362 1350 1353 1359 1364 ...
> ..@ dims : Named int [1:18] 1 500 2 100 1 1 0 0 2 5 ...
> .. ..- attr(*, "names")= chr [1:18] "nt" "n" "p" "q" ...
> ..@ fixef : num [1:2, 1:1000] 1 1 1.093 0.989 1.038 ...
> .. ..- attr(*, "dimnames")=List of 2
> .. .. ..$ : chr [1:2] "(Intercept)" "time"
> .. .. ..$ : NULL
> ..@ nc : int 1
> ..@ ranef : num[1:100, 0 ]
> ..@ sigma : num [1, 1:1000] 0.977 0.853 0.84 0.862 0.921 ...
>
> ##then I did the following which is not in the help file
> > xyplot(x)##check
>
> > x <- mcmcsamp(obj = m41, n = 1e3)
>
> > summary(t(x at fixef))
> (Intercept) emosympt schprob totalnetscr
> Min. : 1.82 Min. :-0.4772 Min. :-0.206 Min. :-0.529
> 1st Qu.:22.11 1st Qu.:-0.1980 1st Qu.: 0.176 1st Qu.:-0.434
> Median :27.10 Median :-0.1144 Median : 0.268 Median :-0.408
> Mean :27.20 Mean :-0.1184 Mean : 0.262 Mean :-0.407
> 3rd Qu.:32.15 3rd Qu.:-0.0416 3rd Qu.: 0.348 3rd Qu.:-0.381
> Max. :51.48 Max. : 0.2662 Max. : 0.636 Max. :-0.264
>
> > colMeans(t(x at fixef))
> (Intercept) emosympt schprob totalnetscr
> 27.200 -0.118 0.262 -0.407
> ##UCL and LCL
> > colMeans(t(x at fixef)) + 1.96*sqrt(colVars(t(x at fixef)))
> (Intercept) emosympt schprob totalnetscr
> 41.537 0.107 0.526 -0.330
>
> > colMeans(t(x at fixef)) - 1.96*sqrt(colVars(t(x at fixef)))
> (Intercept) emosympt schprob totalnetscr
> 12.86248 -0.34386 -0.00259 -0.48348
>
> Respectfully,
>
> Frank R. Lawrence
>
> # -----Original Message-----
> # From: Andrew Robinson [mailto:A.Robinson at ms.unimelb.edu.au]
> # Sent: Saturday, December 06, 2008 2:45 PM
> # To: Frank Lawrence
> # Subject: Re: [R-sig-ME] mcmcsamp
> #
> # Hi Frank,
> #
> # can you provide a minimal, executable example?
> #
> # Cheers
> #
> # Andrew
> #
> # On Fri, Dec 05, 2008 at 02:53:34PM -0500, Frank Lawrence wrote:
> # > I was attempting to run mcmcsamp on an lmer model without success. From
> # the
> # > archive I noted that some users had a similar difficulty a couple of
> months
> # > ago with obtaining fixed effect estimates. I was wondering if there is
> any
> # > new information on using mcmcsamp to obtain confidence intervals for
> fixed
> # > effects from an lmer object.
> # >
> # > Windows Vista, Home Premium. R-2.8
> # >
> # > Respectfully,
> # >
> # > Frank R. Lawrence
> # >
> # > _______________________________________________
> # > R-sig-mixed-models at r-project.org mailing list
> # > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> #
> # --
> # Andrew Robinson
> # Department of Mathematics and Statistics Tel: +61-3-8344-6410
> # University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
> # http://www.ms.unimelb.edu.au/~andrewpr
> # http://blogs.mbs.edu/fishing-in-the-bay/
--
Andrew Robinson
Department of Mathematics and Statistics Tel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
More information about the R-sig-mixed-models
mailing list