[R-sig-ME] mcmcsamp
Frank Lawrence
frl2 at psu.edu
Mon Dec 8 21:09:28 CET 2008
Hi Andrew:
In the past I have been able to obtain confidence intervals for the
parameter estimates using something like the following:
set.seed(281)
nn <- 2e3
ss1 <- mcmcsamp(obj = x, n = nn, v = F)#markov chain sampling from
posterior distribution of parameter estimates
k <- as.matrix(t(apply(X = ss1, MARGIN = 2, FUN = quantile, p =
c(0.025, 0.5, 0.975), na = T, names = T, type = 7)))
colnames(k) <- c('2.5%', '50%', '97.5%')
The series of commands no longer works. I transitioned to something like I
had illustrated below but was not sure it was the most effective or
efficient syntax. I was wondering if there was a better alternative.
Respectfully,
Frank R. Lawrence
# -----Original Message-----
# From: Andrew Robinson [mailto:A.Robinson at ms.unimelb.edu.au]
# Sent: Monday, December 08, 2008 2:52 PM
# To: Frank Lawrence
# Cc: R-sig-mixed-models at r-project.org
# Subject: Re: [R-sig-ME] mcmcsamp
#
# 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