MCMCglmm for binomial models?

Bob,

Just wanted to post the results from my project.  It is still very
much alpha quality software, but the speed is good.  Your sample model
runs in 22secs on my workstation.  I ran the model in JAGS earlier but
it took ages (20min or so).

These are my results which are closer to the MCMCglmm run than the
WinBUGS run.  Eventually, I'll try to make it easy to run these models
directly from R.

The full code for this model is posted in-line below.  However, the
guts are just two functions which just separate the winbugs model into
separate "update" and "logp" functions:

  void update() {
    phi.value = b0.value + b_period2.value*period2 +
b_period3.value*period3 + b_period4.value*period4 +
sum(permutation_matrix*b_herd.value,1) + overdisp.value;
    phi.value = 1/(1+exp(-phi.value));
    sigma_overdisp.value = 1/sqrt(tau_overdisp.value);
    sigma_b_herd.value = 1/sqrt(tau_b_herd.value);
  }
  double logp() const {
    return b0.logp() + b_period2.logp() + b_period3.logp() +
b_period4.logp() + tau_overdisp.logp() + tau_b_herd.logp() +
b_herd.logp(0, tau_b_herd.value) +
      overdisp.logp(0,tau_overdisp.value) + likelihood.logp(size,phi.value);
  }

The speed of the model could probably be improved by consolidating all
the b's into a vector.  I'll try that later on.  Code for this model,
and some other examples will go up on github sometime in the near
future.

I've also cc'd Chris Fonnesbeck, since he may want to chime in on the
pymc front.  I haven't had time to code up an equivalent pymc model,
but it should be very easy to do.

Feedback welcome.

-Whit

results:
warmstrong at krypton:~/dvl/c++/CppBugs/test$ time ./a.out
samples: 17999
b0: -1.52515
b_period2: -1.21971
b_period3: -1.32685
b_period4: -1.88475
tau_overdisp: 1.59118
tau_b_herd: 25.9476
sigma_overdisp: 0.886191
sigma_b_herd: 0.258111
b_herd:
   0.1623
  -0.0524
   0.0852
   0.0080
  -0.0506
  -0.0953
   0.1859
   0.0835
  -0.0690
  -0.0897
   0.0096
  -0.0183
  -0.1496
   0.1133
  -0.1061

real	0m21.930s
user	0m21.920s
sys	0m0.000s
warmstrong at krypton:~/dvl/c++/CppBugs/test$

model code:
class HerdModel: public MCModel {
  const ivec incidence;
  const ivec size;
  const ivec herd;
  const vec period2;
  const vec period3;
  const vec period4;
  int N, N_herd;
  mat permutation_matrix;

public:
  NormalStatic<double> b0;
  NormalStatic<double> b_period2;
  NormalStatic<double> b_period3;
  NormalStatic<double> b_period4;
  UniformStatic<double> tau_overdisp;
  UniformStatic<double> tau_b_herd;
  Deterministic<double> sigma_overdisp;
  Deterministic<double> sigma_b_herd;
  Normal<vec> b_herd;
  Normal<vec> overdisp;
  Deterministic<vec> phi;
  Binomial<ivec> likelihood;

  HerdModel(const ivec& incidence_,const ivec& size_,const ivec& herd_,
            const vec& period2_,const vec& period3_,const vec&
period4_, int N_, int N_herd_):
    incidence(incidence_),size(size_),herd(herd_),
    period2(period2_),period3(period3_),period4(period4_),
    N(N_),N_herd(N_herd_),permutation_matrix(N,N_herd),
    b0(0,0,0.001),b_period2(0,0,0.001),b_period3(0,0,0.001),b_period4(0,0,0.001),
    tau_overdisp(1,0,1000),tau_b_herd(1,0,100),
    sigma_overdisp(1),sigma_b_herd(1),
    b_herd(randn<vec>(N_herd_)),overdisp(randn<vec>(N)),
    phi(randu<vec>(N)),
    likelihood(incidence_,true)
  {
    permutation_matrix.fill(0.0);
    for(uint i = 0; i < herd.n_elem; i++) {
      permutation_matrix(i,herd[i]) = 1.0;
    }
    add(b0);
    add(b_period2);
    add(b_period3);
    add(b_period4);
    add(tau_overdisp);
    add(tau_b_herd);
    add(sigma_overdisp);
    add(sigma_b_herd);
    add(b_herd);
    add(overdisp);
    add(phi);
    add(likelihood);
  }

  void update() {
    phi.value = b0.value + b_period2.value*period2 +
b_period3.value*period3 + b_period4.value*period4 +
sum(permutation_matrix*b_herd.value,1) + overdisp.value;
    phi.value = 1/(1+exp(-phi.value));
    sigma_overdisp.value = 1/sqrt(tau_overdisp.value);
    sigma_b_herd.value = 1/sqrt(tau_b_herd.value);
  }
  double logp() const {
    return b0.logp() + b_period2.logp() + b_period3.logp() +
b_period4.logp() + tau_overdisp.logp() + tau_b_herd.logp() +
b_herd.logp(0, tau_b_herd.value) +
      overdisp.logp(0,tau_overdisp.value) + likelihood.logp(size,phi.value);
  }
};

MCMCglmm for binomial models?

Thread (5 messages)