R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Hierarchical Normal Means Regression Model

hierMeanReg

R Documentation

Hierarchical Normal Means Regression Model

Description

fits a hierarchical normal model of the form E[y_{ij}] = μ_{j} + β_{1}x_{i1}+…+β_{p}x_{ip}

Usage

hierMeanReg(design, priorTau, priorPsi, priorVar,
            priorBeta = NULL, steps = 1000, startValue = NULL,
            randomSeed = NULL)

Arguments

`design`	a list with elements y = response vector, group = grouping vector, x = matrix of covariates or NULL if there are no covariates
`priorTau`	a list with elements tau0 and v0
`priorPsi`	a list with elements psi0 and eta0
`priorVar`	a list with elements s0 and kappa0
`priorBeta`	a list with elements b0 and bMat or NULL if x is NULL
`steps`	the number of Gibbs sampling steps to take
`startValue`	a list with possible elements tau, psi, mu, sigmasq and beta. tau, psi and sigmasq must all be scalars. mu and beta must be vectors with as many elements as there are groups and covariates respectively
`randomSeed`	a random seed for the random number generator

Value

A data frame with variables:

`tau`	Samples from the posterior distribution of tau
`psi`	Samples from the posterior distribution of psi
`mu`	Samples from the posterior distribution of mu
`beta`	Samples from the posterior distribution of beta if there are any covariates
`sigmaSq`	Samples from the posterior distribution of σ^2
`sigma`	Samples from the posterior distribution of sigma

Examples

priorTau <- list(tau0 = 0, v0 = 1000)
priorPsi <- list(psi0 = 500, eta0 = 1)
priorVar <- list(s0 = 500, kappa0 = 1)
priorBeta <- list(b0 = c(0,0), bMat = matrix(c(1000,100,100,1000), nc = 2))

data(hiermeanRegTest.df)
data.df <- hiermeanRegTest.df
design <- list(y = data.df$y, group = data.df$group,
               x = as.matrix(data.df[,3:4]))
r<-hierMeanReg(design, priorTau, priorPsi, priorVar, priorBeta)

oldPar <- par(mfrow = c(3,3))
plot(density(r$tau))
plot(density(r$psi))
plot(density(r$mu.1))
plot(density(r$mu.2))
plot(density(r$mu.3))
plot(density(r$beta.1))
plot(density(r$beta.2))
plot(density(r$sigmaSq))
par(oldPar)

## example with no covariates
priorTau <- list(tau0 = 0, v0 = 1000)
priorPsi <- list(psi0 = 500, eta0 = 1)
priorVar <- list(s0 = 500, kappa0 = 1)

data(hiermeanRegTest.df)
data.df <- hiermeanRegTest.df
design <- list(y = data.df$y, group = data.df$group, x = NULL)
r<-hierMeanReg(design, priorTau, priorPsi, priorVar)

oldPar <- par(mfrow = c(3,2))
plot(density(r$tau))
plot(density(r$psi))
plot(density(r$mu.1))
plot(density(r$mu.2))
plot(density(r$mu.3))
plot(density(r$sigmaSq))
par(oldPar)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Bolstad2)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Bolstad2/hierMeanReg.Rd_%03d_medium.png", width=480, height=480)
> ### Name: hierMeanReg
> ### Title: Hierarchical Normal Means Regression Model
> ### Aliases: hierMeanReg
> 
> ### ** Examples
> 
> priorTau <- list(tau0 = 0, v0 = 1000)
> priorPsi <- list(psi0 = 500, eta0 = 1)
> priorVar <- list(s0 = 500, kappa0 = 1)
> priorBeta <- list(b0 = c(0,0), bMat = matrix(c(1000,100,100,1000), nc = 2))
> 
> data(hiermeanRegTest.df)
> data.df <- hiermeanRegTest.df
> design <- list(y = data.df$y, group = data.df$group,
+                x = as.matrix(data.df[,3:4]))
> r<-hierMeanReg(design, priorTau, priorPsi, priorVar, priorBeta)
           N         mean        stdev        sterr        min           q1
tau     1000    16.236916 3.427045e+01    1.0837267 -89.105729    -6.148252
psi     1000 95491.532177 1.094977e+05 3462.6213807   3.147650 35980.402902
mu.1    1000   276.205849 2.704927e+01    0.8553731 -20.296101   259.308444
mu.2    1000   280.449394 2.625457e+01    0.8302424 -23.738044   265.316729
mu.3    1000   285.409805 2.684243e+01    0.8488322 -23.274001   270.018237
beta.1  1000     5.175377 3.370099e+00    0.1065719  -4.826466     2.930373
beta.2  1000     3.200292 6.782633e+00    0.2144857 -22.029793    -1.556558
sigmaSq 1000  2331.845445 6.968921e+02   22.0376647   3.664022  1830.653794
sigma   1000    47.779282 7.002474e+00    0.2214377   1.914164    42.786140
                 med           q3          max
tau        15.503716 3.843687e+01 2.430817e+02
psi     60723.342950 1.115987e+05 1.780386e+06
mu.1      277.140435 2.928523e+02 3.640613e+02
mu.2      280.813180 2.963322e+02 3.780218e+02
mu.3      286.253053 3.009222e+02 3.719192e+02
beta.1      5.166732 7.284136e+00 1.568636e+01
beta.2      3.507066 7.567394e+00 2.580866e+01
sigmaSq  2213.253863 2.696376e+03 6.531113e+03
sigma      47.045232 5.192664e+01 8.081530e+01
> 
> oldPar <- par(mfrow = c(3,3))
> plot(density(r$tau))
> plot(density(r$psi))
> plot(density(r$mu.1))
> plot(density(r$mu.2))
> plot(density(r$mu.3))
> plot(density(r$beta.1))
> plot(density(r$beta.2))
> plot(density(r$sigmaSq))
> par(oldPar)
> 
> ## example with no covariates
> priorTau <- list(tau0 = 0, v0 = 1000)
> priorPsi <- list(psi0 = 500, eta0 = 1)
> priorVar <- list(s0 = 500, kappa0 = 1)
> 
> data(hiermeanRegTest.df)
> data.df <- hiermeanRegTest.df
> design <- list(y = data.df$y, group = data.df$group, x = NULL)
> r<-hierMeanReg(design, priorTau, priorPsi, priorVar)
           N         mean        stdev        sterr         min           q1
tau     1000     14.09161     32.99325    1.0433380 -78.0954533    -7.507234
psi     1000 138543.64296 263269.76249 8325.3208851   0.4633706 50223.670927
mu.1    1000    310.74755     29.10791    0.9204731 -13.1242793   302.269645
mu.2    1000    316.34499     29.81930    0.9429691 -12.7696201   308.207586
mu.3    1000    321.12488     30.07841    0.9511627 -13.5891102   313.141948
sigmaSq 1000   3048.35491   5410.18680  171.0851286  13.3638613  1989.612778
sigma   1000     51.41741     20.12484    0.6364035   3.6556615    44.605074
                med           q3          max
tau        12.54687     34.42228     118.4996
psi     86368.98897 141996.00599 6510198.9250
mu.1      313.55827    323.78471     365.6160
mu.2      319.79025    328.88741     379.7092
mu.3      323.70052    334.27070     374.2558
sigmaSq  2362.28523   2856.30724   73760.3392
sigma      48.60335     53.44443     271.5885
> 
> oldPar <- par(mfrow = c(3,2))
> plot(density(r$tau))
> plot(density(r$psi))
> plot(density(r$mu.1))
> plot(density(r$mu.2))
> plot(density(r$mu.3))
> plot(density(r$sigmaSq))
> par(oldPar)
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>