[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