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);
?}
};



On Mon, Sep 20, 2010 at 5:47 PM, Bob Farmer <farmerb at gmail.com> wrote:

Thanks for the further detail. ?I wasn't familiar with these two
methods of overdispersion; I think I only understood them in the
WinBUGS sense, which, I believe, is additive (e.g. overdisp[i] as a
parameter). ?I'll have to explore the family() helpfiles a bit more.

Below is a summary of the model outputs from the three methods
described earlier. ?Note that here, I bumped the WinBUGS iterations to
750E3 to reduce the Rhats a bit further.
I also use an adaptation of the bugsParallel function to make my
WinBUGS runs 3x faster (time for this run: ?5.8 minutes); I can
provide details if anybody else would like to use this function (and
some other bugs summary functions I've written).
Hope this is useful or interesting -- I appreciate this discussion!

--Bob

**glmer**

summary(gm3)@coefs

? ? ? ? ? ? ?Estimate Std. Error ? t value
(Intercept) -1.3985351 0.01698321 -82.34810
period2 ? ? -0.9923347 0.02275838 -43.60304
period3 ? ? -1.1286754 0.02429833 -46.45075
period4 ? ? -1.5803739 0.03195596 -49.45474

**MCMCglmm**

summary(mc3)

?Iterations = 749351
?Thinning interval ?= 100001
?Sample size ?= 1000
?DIC: 539.7889
?G-structure: ?~herd
? ? post.mean ?l-95% CI u-95% CI eff.samp
herd ? ?0.2894 1.309e-06 ? ?1.035 ? ? 1000
?R-structure: ?~units
? ? ?post.mean l-95% CI u-95% CI eff.samp
units ? ? ?0.93 0.003825 ? ? 1.98 ? ? 1000
?Location effects: cbind(incidence, size - incidence) ~ period
? ? ? ? ? ?post.mean l-95% CI u-95% CI eff.samp ?pMCMC
(Intercept) ? -1.5263 ?-2.2221 ?-0.9602 ? 1000.0 <0.001 ***
period2 ? ? ? -1.2407 ?-2.2072 ?-0.2449 ? ?947.2 ?0.012 *
period3 ? ? ? -1.3574 ?-2.4534 ?-0.3472 ? 1000.0 ?0.012 *
period4 ? ? ? -1.9126 ?-3.2662 ?-0.8298 ? 1136.0 ?0.002 **

**WinBUGS**

bug3$summary[which(rownames(bug3$summary) %in% params), c("mean", "sd", "2.5%", "97.5%", "Rhat", "n.eff")]

? ? ? ? ? ? ? ? ? ? mean ? ? ? ?sd ? ? ? ?2.5% ? ? ?97.5% ? ? Rhat n.eff
B.0 ? ? ? ? ? ?-1.6070301 0.3556674 -2.33980000 -0.9947425 1.039548 ? ?69
B.period2 ? ? ?-1.2761798 0.4895902 -2.14292500 -0.2391200 1.018612 ? 210
B.period3 ? ? ?-1.5100405 0.6337780 -2.63680000 -0.2266425 1.123427 ? ?22
B.period4 ? ? ?-2.0678331 0.6539004 -3.31892500 -0.8586700 1.018622 ? 270
sigma.overdisp ?1.1975969 0.3147531 ?0.69358404 ?1.8740000 1.045062 ? ?55
sigma.b.herd ? ?0.5120043 0.3825848 ?0.01658429 ?1.3659249 1.542432 ? ? 7
(note that I can provide a link to parameter histograms, if you want
more detail)

--Bob