Laurie Spiegel’s Algorithm

In her article Sonic Set Theory: A Tonal Music Theory for Computers, Laurie Spiegel described a little algorithm, a chord sequence generator, to demonstrate principles of algorithmic composition. Curious of if it could be done in R and how it would sound I wrote a simple sine- and envelope generator and implemented the algorithm in R. Here it is:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
library("audio")
 
# gen_sin(): generate sine signal
#
# sample rate=44100 Hz
# f:    frquency in Hz
# ms:   duration in milliseconds
gen_sin<-function(f, ms)
{
  x<-NULL
  k=44100/(2*pi*f)
  x<-sin(0:(44100*(ms/1000))/k)
  return(x)
}
 
# Check
gen_sin(440,10)
 
# envelope(): generate simple linear AD envelope
#
# x:            signal returned by gen_sin
# attack_ms:    attack time in milliseconds
# decay:        decay time in milliseconds
# max_vol:      maximum volume should not exceed 1 
#               to avoid clipping
envelope<-function(x, attack_ms, decay_ms, max_vol)
{
  length=length(x)
 
  attack_start=1
  attack_end=44100*(attack_ms/1000)
  attack_slope=max_vol/(attack_end-attack_start)
 
  decay_end=length(x)
  decay_start=decay_end-44100*(decay_ms/1000)
  decay_slope=max_vol/(decay_end-decay_start)
 
  env<-rep(max_vol, length)
 
  volume_attack=0
  for(i in attack_start:attack_end)
  { 
    env[i]=volume_attack
    volume_attack=volume_attack+attack_slope
  }
  volume_decay=max_vol
  for(i in (decay_start):decay_end)
  {
    env[i]=volume_decay
    volume_decay=volume_decay-decay_slope
  }
  return(env)  
}
 
# Checks
x<-gen_sin(440,100)
length(x)
plot(x)
env<-envelope(x,50,50,0.5)
length(env)
plot(env)
 
# gen note(): generate note
#
# f:          frequency in Hertz
# ms:         duration in milliseconds
# attack_ms:    attack time in milliseconds
# decay:        decay time in milliseconds
# max_vol:      maximum volume should not exceed 1 
#               to avoid clipping
gen_note<-function(f, ms, attack_ms, decay_ms, max_vol)
{
  x<-gen_sin(f,ms)
  env<-envelope(x,attack_ms, decay_ms, max_vol)
  env_signal<-x*env
  return(env_signal)
}
 
x<-gen_note(440,1000,100,500,1)
play(x)
plot(x)
 
# gen_triad(): generate frequencies of a triad
# note:        base note of triad
# freq_scale:  scale frequencies
gen_triad<-function(note, freq_scale)
{
  chord<-NULL
  for(i in 1:3)
  {
    index=((note-1)%%16)+1 # R starts counting with 1
    note=note+2
    chord<-c(chord, gen_note(freq_scale[index], 400,40,200,0.8))
  }
    return(chord)
}
 
# Frequencies of the C major scale
freq_scale<-c(261.63, 293.66,  # C4 D4 
              329.63, 349.23,  # E4 F4
              392.00, 440.00,  # G4 A4
              493.88, 523.25,  # B4 C5
              587.33, 659.25,  # D5 E5
              698.46, 783.99)  # F5 G5  
 
# Check
plot(freq_scale)
# Play C Major chord
x<-gen_triad(1, freq_scale)
play(x)
 
### Generate piece ###
# Source of the alogrithm:
# Spiegel, L. (1982). Sonic Set Theory: A Tonal Music Theory for Computers. 
# In Proceedings of the Second Annual Symposium on Small Computers and the Arts.
#     
stages<-matrix(c(0, 0, 1, 0, 0, 0, 0,
                 0.5, 0, 0, 0, 0, 0.5, 0,
                 0, 0.5, 0, 0.5, 0, 0, 0,
                 0, 0, 0, 0, 0.5, 0, 0.5,
                 0.5, 0, 0, 0, 0, 0.5, 0), byrow=TRUE, ncol=7)
 
# Base notes of chords
base_notes<-rep(1:7) 
 
piece<-NULL
cycles=0
start=1
max_cycles=4
 
while(cycles<max_cycles)
{
  for(i in start:5)
  {
    index<-sample(base_notes,1, p=stages[i,])    
    piece<-c(piece, gen_triad(index, freq_scale))
    cat(index)
  }
  start=sample(2:5,1)
  cycles=cycles+1
}
 
play(piece)

Hope I got it right. Read the article, play around with the probabilities and have fun! You do not have to use C major of course. Oh, and be careful with your speakers.

ProTrackR

While exploring a bit if R could be used for algorithmic composition I stumbled over the R ProTrackR package.

Check out this funky tune:

1
2
3
require("ProTrackR")
x<-read.module("https://api.modarchive.org/downloads.php?moduleid=40475#ELYSIUM.MOD")
playMod(x)

By the way, algorithmic composition in R? I think so, thanks to the R audio package.

Wow, I have a blog?

Almost three years have passed since I last updated this blog, a long time!
What happened in the meantime? Well, I have been working a bit in the domain of aesthetics and thought a little bit about Egon Brunswik’s probabilistic functionalism among other things.

Tomorrow, I am leaving for Vienna to get some feedback with regards to an R-package.

Well, feedback from friends with regard to this blog was more like “I don’t understand a word!” or “Can’t you just post pictures???”.

Ok. Photography got neglected by me for quite some time, but here are some more recent shots I like.

This is a shot of Marzahn. Now, everyone seems to have an opinion about Marzahn. If you get a chance, try it out by yourself.


This is Stefan, a good friend. He is a sysadmin. Apart from that, he has a good eye for forms and loves taking high-contrast black and white macro shots of plants on film.


The botanical garden of Berlin is a good place if you are interested in plants. Sometimes it makes me feel like being in a space station.


And this is me in front of a painting by Hieronymus Bosch taken by Lem.

JAGS as a tool for simulation studies

In one of my last posts I hinted to the fact that JAGS could be used for simulation studies. Let’s have a look! So basically simulation studies are very useful if one wants to check if a model one is developing makes sense. The main aim is to assess, if the parameters one has plugged into a data simulation process based on the model could be reasonably recovered by an estimation process. In JAGS it is possible to simulate data given a model and a set of known parameters and to estimate the parameters given the model and the simulated data set. For today’s example I have chosen a simple psychometric latent variable measurement model from classical test theory, which is very well understood and can be easily estimated with most software capable of latent variable modeling. Here it is:

\mu_{1i}=\lambda_1\cdot\eta_i+\beta_1
\mu_{2i}=\lambda_2\cdot\eta_i+\beta_2
\mu_{3i}= 1\cdot\eta_i+\beta_3.

y_{1i} \sim \mbox{dnorm} (\mu_{1i}, \sigma^2_{e1})
y_{2i} \sim \mbox{dnorm} (\mu_{2i}, \sigma^2_{e2})
y_{3i} \sim \mbox{dnorm} (\mu_{3i}, \sigma^2_{e3})

\eta_i \sim \mbox{dnorm} (0, \sigma^2_{\eta})

The basic idea is that a latent trait \eta is measured by three items (y_{1i} to y_{3i}) with different difficulties (\beta_1 to \beta_3), different loadings (\lambda_1 and \lambda_2) and different measurement error variances (\sigma^2_{e1} to \sigma^2_{e3}). For identification reasons, the loading of item 3 is set to 1. Furher the variance of the trait (\sigma^2_{\eta}) has to be estimated from data.

For simulating data with JAGs we have to specify a data block. The model specification is straight forward, note that we have to hand over the true parameters to JAGS, for example with an R script. This is pretty handy because one could design a thorough simulation study for different sample size and parameter combinations. The estimation part has to be specfied in a model block.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
#Simulation
data{
 
  for(i in 1:n)
  {
      y1[i]~dnorm(mu1.true[i], tau.y1.true);
      y2[i]~dnorm(mu2.true[i], tau.y2.true);
      y3[i]~dnorm(mu3.true[i], tau.y3.true);
 
      mu1.true[i]<-lambda1.true*eta.true[i]+beta1.true
      mu2.true[i]<-lambda2.true*eta.true[i]+beta2.true
      mu3.true[i]<-eta.true[i]+beta3.true
  }
 
  tau.y1.true<-pow(sigma.y1.true, -2); 
  tau.y2.true<-pow(sigma.y2.true, -2);
  tau.y3.true<-pow(sigma.y3.true, -2);
 
}
 
# Estimation
model{
 
  for(i in 1:n)
  {
    y1[i]~dnorm(mu1[i], tau.y1);
    y2[i]~dnorm(mu2[i], tau.y2);
    y3[i]~dnorm(mu3[i], tau.y3);
 
    mu1[i]<-lambda1*eta[i]+beta1
    mu2[i]<-lambda2*eta[i]+beta2
    mu3[i]<-eta[i]+beta3
  }
 
  for(i in 1:n)
  {
    eta[i]~dnorm(0,tau.eta);
  }
 
  # Priors for the loadings
  lambda1~dnorm(0,1.0E-3)
  lambda2~dnorm(0,1.0E-3)
  bias_lambda1<-lambda1-lambda1.true
 
  # Priors for the easiness parameters
  beta1~dnorm(0,1.0E-3)
  beta2~dnorm(0,1.0E-3)
  beta3~dnorm(0,1.0E-3)
  bias_beta1<-beta1-beta1.true
 
  # Hyperprior for latent trait distribution
  tau.eta<-pow(sigma.eta,-2)
  sigma.eta~dunif(0,10) 
 
  # Priors for the measurement errors
   tau.y1<-pow(sigma.y1, -2); 
   tau.y2<-pow(sigma.y2, -2);
   tau.y3<-pow(sigma.y3, -2);
   sigma.y1~dunif(0,10);
   sigma.y2~dunif(0,10); 
   sigma.y3~dunif(0,10); 
}

This may look a bit daunting, but it really took only 15 minutes to set that up.

The simulation is controlled by an R script, which hands over the true parameters to JAGS:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
library(R2jags)
library(mcmcplots)
 
n=500 # number of subjects
 
# generate 'true' latent trait values
eta.true<-rnorm(n, 0, 4)
# generate 'true' factor loadings
lambda1.true=1
lambda2.true=0.5
 
beta1.true=1
beta2.true=2
beta3.true=3
 
# generate 'true' error standard deviations
sigma.y1.true=1
sigma.y2.true=0.5
sigma.y3.true=0.75
 
data<-list("n", "eta.true", "lambda1.true", "lambda2.true", "sigma.y1.true", "sigma.y2.true", "sigma.y3.true",
           "beta1.true", "beta2.true", "beta3.true")
 
parameters<-c("lambda1","lambda2", "sigma.eta", "sigma.y1", "sigma.y2", "sigma.y3", 
              "beta1", "beta2", "beta3", "bias_beta1")
 
output<-jags.parallel(data, inits= NULL, parameters, model.file="simple_1f.txt",n.chains=2, n.iter=35000, n.burnin=30000)
attach.jags(output)
denplot(output)
traplot(output)
plot(output)

If we inspect the posterior densities based on one dataset with sample size 500 we see that we pretty much hit the spot. That is, all posterior densities are covering the generating true parameters.
JAGS_sim
Please note that in the estimation code a posterior density for the difference between the data generating parameter \beta_1 and its’ samples from the posterior is specified. The zero is very well covered, but keep in mind that this is based only on one dataset. To replicate a more traditional approach to bias, a viable route could be to simulate several hundreds of datasets and to record for example the posterior mean for each parameter and data set to assess bias.

As always, keep in mind that this is pretty experimental code and estimating latent variable models with JAGS is not very well explored yet as far as I know. However, the main point of this post was to show how JAGS could be used in combination with R to run simulation studies in a probabilistic graphical framework. Should work for multilevel models or generalized mixed models too.

Edit: I just stumbled about this recent blog post by Robert Grand. So stan should work too. Interesting!

Pictures!

This blog looks super technical right now and one may get the impression that I am only munching probability distributions day in and day out! No! I am also munching pizza!

Keep in mind that this is my personal blog so I may also post things that interest me apart from statistics. One of those things are pictures!

There are quite some photos by me online but not so many in which I am depicted. So here is a recent one of me making pizza dough at a friends’ place:

georg_making_pizza _dough

Picture by Martin.

 

Bayesian estimation with JAGS and MPlus: a quick practical comparison

Bayesian estimation techniques were added to Mplus version 7 some time ago. Of course I was interested if MPlus and JAGS come to comparable results when let loose on data. 😉

In particular I was interested in the crossclassified module of MPlus, so I wrote some code to simulate data with R and subsequently recovered the parameters with JAGS and MPlus. Recently I found out how JAGS could also be used as a data simulator for simulation studies, but this is a different topic.

So graphically, the model I used for simulation looks like this:
Diagram1
I have omitted the level-1 residual variances and the variances of the latent distributions in the figure for simplicity reasons. The full model specifications can be gleaned from the source code attached to this post.

To make a long story shot: both approaches, JAGs and MPlus, pleasently come to very comparable results with regard to the parameters. At the moment I am too busy with other things to harp on this more, but I thought it is useful to put the code on the blog in case someone is interested to work on a similar topic in more detail. Thanks to Mr Asparouhov for clarifying some questions with regard to MPlus.

The files:
test.R: code to simulate data from the model, recover the parameters with JAGs, also a bit of convergence diagnostics

test.txt: JAGs model file

cc.inp: MPlus input file

cc.dat: Simulated data. Will be overwritten once the simulation part of test.R is executed.

Get the files here. The code is not very polished but should give an idea how to specify the model in MPlus and JAGs.

Blogroll

The blog has been neglected for quite some time. At the moment I am using it more or less as a public notepad, but comments are definately invited! However, I added a blogroll.

The book by Gelman and Hill on multilevel hierarchical models helped me quite a bit with regard to the topic of multilevel/hierarchical models and it is still very useful, especially the part about modeling multivariate normal distributions for random effects with MCMC.

Felix maintains quite an active and interesting blog with topics related to science, statistical methods and R.

John R. Kruschke is known for his book on doing Bayesian data analysis, as far as I know one of the first in the field for psychologists. A book review by Wolf Vanpaemel and Francis Tuelinckx appeared in the Journal of Mathematical Psychology.

R-Bloggers is a central hub of content collected from bloggers who write about R.

Rotational Invariance in factor models II

I have been working a bit more on estimating simple factor models with JAGS. I posted the problem from the last post over there at the open discussion forum and had some very helpful input from Vincent Arel-Bundock from Michigan University and one of his co-authors, Walter R. Mebane which ultimately led to the code of the current post.

Here I wanted to know how to estimate a simple factor model for continuous manifest variables and three latent, intercorrelated factors. The code is simple and straight forward. I put a scaled inverse Wishart prior on the variance-covariance matrix of factor loadings, following the approach in Gelman & Hill (2007, p.376). For solving the invariance problem I set the loadings for one item on each factor to 1. In addition, the means of the latent trait distributions (factor values) were all constrained to zero. This is a very traditional approach, that is also implemented in e.g. MPlus, the STATA program GLLAMM or the R package lavaan.

In the following code I use the Holzinger & Swineford dataset from the lavaan tutorial for comparison. First I estimate the model with lavaan and in a second step with JAGS. Both approaches show very comparable results, but keep in mind that this is somewhat experimental. When comparing coefficients please note that I coded the model in such a way that JAGS returns estimated standard deviations and not variances for the manifest items and latent variables.

Ok, enough talk, here is the code. First the R-script:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
## holzinger.R
# Experimental!
 
library(lavaan)
 
# From the lavaan tutorial
HS.model <- ' visual  =~ x1 + x2 + x3      
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '
 
fit <- cfa(HS.model, data=HolzingerSwineford1939, meanstructure=TRUE)
 
summary(fit, standardized=TRUE)
 
library(R2jags)
data<-HolzingerSwineford1939
## Transform data to long format----
#  the complicated way
 
# Dependent variable
y<-c(data$x1, data$x2, data$x3, data$x4, data$x5, data$x6, data$x7, data$x8, data$x9)
 
# Item identifier
item<-c(rep(1, 301), rep(2, 301), rep(3, 301), rep(4, 301), rep(5, 301), rep(6, 301), rep(7, 301), rep(8, 301), rep(9, 301))
 
# Construct identifier (1: visual, 2: textual, 3: speed)
const<-c(rep(1,3*301 ), rep(2,3*301 ), rep(3,3*301 ))
 
# Person id
id<-rep(seq(1:301),9)
## ---------------------------------
# Fitting the model with JAGS
W=(diag(3))
 
dat<-list("y", "item", "const", "id", "W")
parameters<-(c("beta", "sigma.y", "lam", "sigma.eta1", "sigma.eta2", "sigma.eta3", "rho_eta13", "rho_eta12", "rho_eta23", "eta"))
 
output<-jags.parallel(dat, inits=NULL, parameters,  model.file="3fa.txt",n.chains=2, n.iter=3000, n.burnin=2000)
plot(output)

Next the model code for JAGS:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
## 3fa.txt
# Experimental!
model{
 
  # Likelihood
  for(i in 1:2709) # Loop over observations
  {
      # Linear decomposition could probably done more
      # elegantly with matrix algebra
      mu[i]<-beta[item[i]] +
         lam[const[i], item[i]]*eta[const[i], id[i]]+
         lam[const[i], item[i]]*eta[const[i], id[i]]+
         lam[const[i], item[i]]*eta[const[i], id[i]]
 
      # Distributional assumptions for the 
      # observed variables
         y[i]~dnorm(mu[i], tau.y[item[i]])
  }
 
  # Priors for the items
    for(i in 1:9)
    {
      beta[i]~dnorm(0,1.0E-3) 
      tau.y[i] <- pow(sigma.y[i], -2)  
      sigma.y[i] ~ dunif (0, 100)   
    }
 
    # Matrix of loadings [construct, item]
    # Constraints and priors 
     lam[1,1]<-1
     lam[2,1]<-0
     lam[3,1]<-0
     lam[1,2]~dnorm(0,1.0E-3) 
     lam[2,2]<-0
     lam[3,2]<-0
     lam[1,3]~dnorm(0,1.0E-3) 
     lam[2,3]<-0
     lam[3,3]<-0
     lam[1,4]<-0
     lam[2,4]<-1
     lam[3,4]<-0
     lam[1,5]<-0
     lam[2,5]~dnorm(0,1.0E-3)
     lam[3,5]<-0
     lam[1,6]<-0
     lam[2,6]~dnorm(0,1.0E-3)
     lam[3,6]<-0
     lam[1,7]<-0
     lam[2,7]<-0
     lam[3,7]<-1
     lam[1,8]<-0
     lam[2,8]<-0
     lam[3,8]~dnorm(0,1.0E-3)
     lam[1,9]<-0
     lam[2,9]<-0
     lam[3,9]~dnorm(0,1.0E-3)
 
    # Scaled inverse Wishart prior for
    # variance-covariance matrix of factor values
    # after Gelman and Hill (2007)
 
    for(i in 1:301)
    {  
      eta[1,i]<-xi.eta1*eta.raw[i,1]
      eta[2,i]<-xi.eta2*eta.raw[i,2]
      eta[3,i]<-xi.eta3*eta.raw[i,3]
      eta.raw[i,1:3]~dmnorm(eta.raw.hat[i,],tau.eta.raw[,])
      eta.raw.hat[i,1]<-mu.eta1.raw
      eta.raw.hat[i,2]<-mu.eta2.raw
      eta.raw.hat[i,3]<-mu.eta3.raw
    }
      mu.eta1<-xi.eta1*mu.eta1.raw
      mu.eta2<-xi.eta2*mu.eta2.raw
      mu.eta3<-xi.eta3*mu.eta3.raw
      mu.eta1.raw<-0
      mu.eta2.raw<-0
      mu.eta3.raw<-0
      xi.eta1~dunif(0,100)
      xi.eta2~dunif(0,100)
      xi.eta3~dunif(0,100)
      tau.eta.raw[1:3,1:3] ~ dwish(W[,],df)
      df<-4
      sigma.eta.raw[1:3,1:3]<-inverse(tau.eta.raw[,])
      sigma.eta1<-xi.eta1*sqrt(sigma.eta.raw[1,1])
      sigma.eta2<-xi.eta2*sqrt(sigma.eta.raw[2,2])
      sigma.eta3<-xi.eta3*sqrt(sigma.eta.raw[3,3])
      rho_eta12<-sigma.eta.raw[1,2]/sqrt(sigma.eta.raw[1,1]*sigma.eta.raw[2,2])
      rho_eta13<-sigma.eta.raw[1,3]/sqrt(sigma.eta.raw[1,1]*sigma.eta.raw[3,3])
      rho_eta23<-sigma.eta.raw[2,3]/sqrt(sigma.eta.raw[2,2]*sigma.eta.raw[3,3])
}

Rotational invariance in factor models

Long time no post. Meanwhile I got my Ph.D. and I am still MCMC’ing. In particular I am interested in estimating latent variable models via MCMC with possible extensions to multilevel latent variable models. To get a grip on it I started with defining a simple one-factorial latent variable model and stumbled over the phenomenon of rotational invariance which is described in quite some detail in this techical report by Erosheva and McKay Curtis. In short, the likelihood surface of factor models is multimodal which makes it difficult to estimate this type of model with Monte Carlo methods. If one is using multiple chains they may operate on different modes of the likelihood surface. This phenomenon is also known as pole switching or label switching in latent class analysis or mixture modeling. Erosheva and McKay Curtis offer a solution based on postprocessing the markov chains. I am wondering if this is feasible, if the models become more complicated.

If you want to have a more hands on experience with regard to the problem, have a look at this JAGS model.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
model{
  # Likelihood
  # i: Persons
  # k: Items
 
  for(i in 1:n)
  {
    for(j in 1:k)
    {
      mu[i,j]<-beta[j]+lambda[j]*eta[i]
      y[i,j]~dnorm(mu[i,j], tau.y[j])
    }
  }
 
  # Priors
  for(j in 1:k)
  {
  # Item-Residuals
    tau.y[j] <- pow(sigma.y[j], -2)
    sigma.y[j] ~ dunif (0, 100)
  # Intercepts
    beta[j]~ dnorm(0,1.0E-3) 
  # Loadings 
   lambda[j]~dnorm(0,1.0E-3)
  }
 
 
  # Factor-Values
  for(i in 1:n)
  {
  # Variance of latent score distribution fixed to 1
    eta[i]~dnorm(0,1)
  }
 
}

In Psychometrics this model is also known as tau congeneric measurement model for continuous observed variables. Here I have fixed the variance of the latent trait distribution to 1 for identifiability reasons.

The following R-Script simply generates some data from the model with know parameters and hands the data over to JAGS.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
library(R2jags)
 
# Number of Persons
n=100
# Number of Items
k=10
# Difficulties
beta.gen=rnorm(k, mean=0, sd=1)
# Abilities
a<-(rnorm(n, 0, sd=1))
eta.gen<-a-mean(a)
 
# Error standard deviations
sigma.gen=runif(k)*3
y<-matrix(rep(0,n*k), nrow=n, ncol=k)
 
# Produce data
for(i in 1:n)
{
  for(j in 1:k)
  {
    y[i,j]=rnorm(1,mean=beta.gen[j]+eta.gen[i], sd=sigma.gen[j])
  }
}
 
data<-list("y", "n", "k")
 
parameters<-c("beta", "sigma.y", "lambda", "eta")
 
# Starting values
inits <- function(){
list(beta = rnorm(k), lambda=c(rnorm(k)), eta=c(rnorm(n)), sigma.y=runif(k)*100)
}
 
output<-jags(data, inits=NULL, parameters, model.file="1f_lat_var.txt",n.chains=2, n.iter=5000, n.burnin=2500)
 
attach.jags(output)
plot(output)

The main interest here is to check, if the data generating parameters could be reasonably recovered. If you run the code – and if you are lucky – the two employed markov chains will converge to one solution. But it may well be that the two chains will end up operating in two different modes of the likelihood surface. In this case you will see a reflection of the two chains for the loadings and the factor values around the y-axis in the plots. This does not mean that the model is ‘wrong’ or something like that. It simply shows that there are two equivalent solutions. In practice one could simply choose one chain and be happy with that, but if one wants to extend the model by introducing more latent variables or even multilevel structures, the multimodality could be problematic when using MCMC-methods. So, if anyone has a simple solution for dealing with this phenomenon, let me know! Meanwhile I am having a look at pymc.

Some short thoughts on MCMC methods

While MCMC is certainly a useful method and has the potential to adress several difficulties that may exist with more traditional methods of parameter estimation, in my opinion it is not a magic bullet for several reasons:

1. The estimation procedure can be slow, especially for larger datasets.
2. The use of prior distributions may deter some researchers.
3. Sometimes it is not clear which prior distribution is appropriate.
4. Model specification in code requires a thorough understanding of the model in question from a formal point of view.

On the other hand, the potential advantages are:

1. Unified handling of a wide scope of models within the framework of graphical models.
2. Unified method for checking a model (posterior predictive checks).
3. Possibility of obtaining posterior distributions of derived statistics, such as effect size indices.
4. Model based handling of missing values and possibility of obtaining posterior distributions of missings in the dependent variable.
5. Possibility of rapid model prototyping.

For building confidence it may be useful to apply the MCMC method to the wide range of models used in the social and behavioral sciences, to cross check against the traditional criteria of parameter estimation (bias, consistency, efficiency) by simulation studies and to examine the influence of the prior distributions on the posterior distributions. But I guess this has already been done. I should check the literature.