[R-sig-ME] mcmcsamp

Mon Dec 8 16:24:05 CET 2008

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/