hbmem (version 0.2)

sampleNormR: Function sampleNormR

Description

Samples posterior of mean parameters of the hierarchical linear normal model with the effects a linear function of some other variable.

Usage

sampleNormR(sample, phi,blockD,y,subj, item, lag, I, J, R, nsub, nitem,
    s2mu, s2a, s2b, meta, metb, sigma2, sampLag)

Arguments

sample

Block of linear model parameters from previous iteration.

y

Vector of data

phi

Vector of linear slopes on effects.

blockD

Block of parameters that will serve as the means of random effects

subj

Vector of subject index, starting at zero.

item

Vector of item index, starting at zero.

lag

Vector of lag index, zero-centered.

I

Number of subjects.

J

Number of items.

R

Total number of trials.

nsub

Vector of length (I) containing number of trials per each subject.

nitem

Vector of length (J) containing number of trials per each item.

s2mu

Prior variance on the grand mean mu; usually set to some large number.

s2a

Shape parameter of inverse gamma prior placed on effect variances.

s2b

Rate parameter of inverse gamma prior placed on effect variances. Setting both s2a AND s2b to be small (e.g., .01, .01) makes this an uninformative prior.

meta

Matrix of tuning parameter for metropolis-hastings decorrelating step on mu and alpha. This hould be adjusted so that .2 < b0 < .6.

metb

Tunning parameter for decorrelating step on alpha and beta.

sigma2

Variance of distribution.

sampLag

Logical. Whether or not to sample the lag effect.

Value

• The function returns a list. The first element of the list is the newly sampled block of parameters. The THIRD element contains a vector of 0s and 1s indicating which of the decorrelating steps were accepted.

References

Not published yet.

See Also

hbmem

Examples

library(hbmem)

I=50
J=100
M=500
B=I+J+4
mu=.5
muS2=0
s2a=.2
s2b=.2
s2aS2=0
s2bS2=0

phi=c(.2,.08)
blockD=rep(0,B)
blockD[2:(I+1)]=rnorm(I,0,.5)
blockD[(I+2):(I+J+1)]=rnorm(J,0,.5)

    R = I * J
    alpha = rnorm(I, phi[1]*blockD[2:(I+1)], sqrt(s2a))
    beta =  rnorm(J, phi[2]*blockD[(I+2):(I+J+1)], sqrt(s2b))
    alphaS2 = rnorm(I, 0, sqrt(s2aS2))
    betaS2 = rnorm(J, 0, sqrt(s2bS2))
    subj = rep(0:(I - 1), each = J)
    item = rep(0:(J - 1), I)
    lag = rep(0, R)
    resp = 1:R
    for (r in 1:R) {
        mean = mu + alpha[subj[r] + 1] + beta[item[r] + 1]
        sd = sqrt(exp(muS2 + alphaS2[subj[r] + 1] + betaS2[item[r] + 1]))
        resp[r] = rnorm(1, mean, sd)
    }
    sim=(as.data.frame(cbind(subj, item, lag, resp)))



blockR=matrix(0,M,B)
blockR[1,c(I+J+2,I+J+3)]=c(.1,.1)
met=c(.1,.1)
b0=c(0,0)

for(m in 2:M)
  {
tmp=sampleNormR(blockR[m-1,],phi,blockD,sim$resp,sim$subj,sim$item,sim$lag,I,J,I*J,table(sim$sub),table(sim$item),10,.01,.01,met[1],met[2],1,1)
blockR[m,]=tmp[[1]]
b0=b0+tmp[[3]]
}

est=colMeans(blockR)

par(defpar(2,3))
plot(blockR[,1],t='l')
abline(h=mu,col="blue")
plot(blockR[,I+J+2],t='l')
abline(h=s2a,col="blue")
plot(blockR[,I+J+3],t='l')
abline(h=s2b,col="blue")

plot(est[2:(I+1)]~alpha);abline(0,1,col="blue")
plot(est[(I+2):(I+J+1)]~beta);abline(0,1,col="blue")

#Compare estimates from regular normal ones:

s.block=matrix(0,nrow=M,ncol=B)
met=c(.1,.1);b0=c(0,0)
for(m in 2:M)
{
tmp=sampleNorm(s.block[m-1,],sim$resp,rep(0,length(sim$resp)),sim$subj,sim$item,sim$lag,1,I,J,R,R,table(sim$subj),
table(sim$item),100,.01,.01,met[1],met[2],1,1)
s.block[m,]=tmp[[1]]
b0=b0 + tmp[[2]]
}

est2=colMeans(s.block)

par(defpar(1,2))
plot(est[2:(I+1)]~est2[2:(I+1)]);abline(0,1,col="blue")
plot(est[(I+2):(I+J+1)]~est2[(I+2):(I+J+1)]);abline(0,1,col="blue")