h1

How-to: VIM, GNU R, and LaTeX

October 24, 2009

Recently I switched from Emacs to VIM. Mostly because I’ve felt that the default keybindings are more sane. The first thing I needed to do was to find a plugin for VIM that worked for R. Well here is that plugin. The follow the directions there for installation. (If you need help leave me a comment). VIM will automatically detect that it is in R mode when opening a *.R file.

Two quick tips to get you started:

# To open R
\rf
# To send code to R
\d

Also, if you need a LaTeX plugin for VIM here it is. It also has very sane keybindings. I recommend not using the Debian package as it didn’t work for me. To get LaTeX to compile in PDF mode by default I added the following to ~/.vimrc:


" Set pdf as default target
let g:Tex_DefaultTargetFormat='pdf'

That’s it. And in complete disclosure, I am really using the GTK-front end of VIM known as GVIM.

h1

SEM in WinBUGS

October 24, 2009

For my MITER research, I have been working with a faculty member in the University of Minnesota’s Institute of Child Development. I have been involved on a couple of projects with this group. One project involved longitudinal modeling of work competency (from childhood to adulthood) using structural equation modeling (SEM). This project examined how childhood competence in externalizing behaviors/conduct, academics, and socializing behaviors predict work as an adult. If you’re interested in this work, it has just been accepted and should be in press soon. (Feel free to email me for a copy of this manuscript and I’ll see if I can get a copy).

To do the SEM, we initially tried the sem package in R. This proved to be to basic for us as it was unable to do a multi-group analysis. So we settled on MPlus. Even in MPlus, we were unable to run the models we really wanted to because the number of integration points were too many and MPlus would just error out. We dealt with it the best we could and I think our manuscript is strong and transparent about our computational problems.

After this paper was submitted, I become more interested in examining what sort of alternatives exist. Is there a Bayesian book on SEM? Well apparently there is. In the Bayesian framework, latent constructs are just dealt with as missing data problems, making them more computationally solvable. Additionally, using Bayesian you work directly from your observed data and not from a covariance matrix. Because of this you can create more highly complex, developmentally appropriate models that are infeasible in traditional software. Even Bollen is writing about this.

To help you get started. I am going to paste the code that I have used to replicate our models in WinBUGS. Unfortunately I can not share the data. But if you have specific data you can share then I’d be happy to modify the code as appropriate.


model {
for(i in 1:N){
##########################
# measurement equation model
##########################
for(j in 1:P){
y[i,j]~dnorm(mu[i,j],psi[j])
}
##########################
# alp[i] corresponds to the intercepts
##########################
mu[i,1]<-xi[i,1]+alp[1] ## Ext1
mu[i,2]<-lam[1]*xi[i,1]+alp[2]
mu[i,3]<-lam[2]*xi[i,1]+alp[3]
mu[i,4]<-xi[i,2]+alp[4] ## Soc1
mu[i,5]<-lam[3]*xi[i,2]+alp[5]
mu[i,6]<-lam[4]*xi[i,2]+alp[6]
mu[i,7]<-xi[i,3]+alp[7] ## Aca1
mu[i,8]<-lam[5]*xi[i,3]+alp[8]
mu[i,9]<-lam[6]*xi[i,3]+alp[9]
mu[i,10]<-xi[i,4]+alp[10] ## Ext2
mu[i,11]<-lam[7]*xi[i,4]+alp[11]
mu[i,12]<-lam[8]*xi[i,4]+alp[12]
mu[i,13]<-xi[i,5]+alp[13] ## Soc2
mu[i,14]<-lam[9]*xi[i,5]+alp[14]
mu[i,15]<-lam[10]*xi[i,5]+alp[15]
mu[i,16]<-xi[i,6]+alp[16] ## Aca2
mu[i,17]<-lam[11]*xi[i,6]+alp[17]
mu[i,18]<-lam[12]*xi[i,6]+alp[18]
mu[i,19]<-xi[i,7]+alp[19] ## Work2
mu[i,20]<-lam[13]*xi[i,7]+alp[20]
mu[i,21]<-lam[14]*xi[i,7]+alp[21]
##########################
# structural equation model
##########################
xi[i,1:3]~dmnorm(u[1:3],phi[1:3,1:3]) ## factor loadings are multivariate normal
xi[i,4]~dnorm(ext2[i],psd)
xi[i,5]~dnorm(soc2[i],psd)
xi[i,6]~dnorm(aca2[i],psd)
xi[i,7]~dnorm(work2[i],psd)
##########################
# continuity Paths
##########################
ext2[i]<-gam[1]*xi[i,1]
soc2[i]<-gam[2]*xi[i,2]
aca2[i]<-gam[3]*xi[i,3]
work2[i]<-gam[4]*xi[i,1] + gam[5]*xi[i,1] + gam[6]*xi[i,1]
} #end of i
##########################
# priors on intercepts
##########################
for(j in 1:21){alp[j]~dnorm(0.0, 1.0)}
##########################
# priors on loadings and coefficients
##########################
lam[1]~dnorm(0.8,psi[1])
lam[2]~dnorm(0.8,psi[2])
lam[3]~dnorm(0.8,psi[3])
lam[4]~dnorm(0.8,psi[4])
lam[5]~dnorm(0.8,psi[5])
lam[6]~dnorm(0.8,psi[6])
lam[7]~dnorm(0.8,psi[7])
lam[8]~dnorm(0.8,psi[8])
lam[9]~dnorm(0.8,psi[9])
lam[10]~dnorm(0.8,psi[10])
lam[11]~dnorm(0.8,psi[11])
lam[12]~dnorm(0.8,psi[12])
lam[13]~dnorm(0.8,psi[13])
lam[14]~dnorm(0.8,psi[14])
for(j in 1:6){gam[j]~dnorm(0.5, psd)}
##########################
#priors on precisions
##########################
for(j in 1:P){
psi[j]~dgamma(9.0, 4.0)
sgm[j]<-1/psi[j]
}
psd~dgamma(9.0, 4.0)
sgd<-1/psd
phi[1:3,1:3]~dwish(R[1:3,1:3], 5)
phx[1:3,1:3]<-inverse(phi[1:3,1:3])
} #end of model

I’d be happy to answer any questions about the code. Note this is just the model, not the starting values. I usually let WinBUGS generate them for me. Bad practice I know, but as long as my chains mix I am not worried.

h1

WinBUGS and OpenBUGS on Debian GNU/Linux Lenny

October 24, 2009

A quick how-to for installing WinBUGS and OpenBUGS on Debian GNU/Linux Lenny. I am sure this will work on other distributions. The only issue might be that certain distributions ship different versions of wine . I am use wine 1.0.1, with other versions of wine YMMV.

First, install wine

su -
aptitude install wine

Second, go here and download the WinBUGS and OpenBUGS binaries. It is important to note that OpenBUGS is the development version of WinBUGS and it’s open-source.

I also encourage you to check out JAGS. I have created a 64-bit package of JAGS that is available here for Debian. If you’re a R user and do a lot of multilevel modeling, I recommend MCMCglmm .

Third, assuming you’ve downloaded these files to ~/Downloads, for WinBUGS do the following:

unzip ~/Downloads/winbugs14.zip

And for OpenBUGS do the following:


mkdir -p ~/Downloads/openbugs
mv ~/Downloads/OpenBUGS.zip ~/Downloads/openbugs
unzip ~/Downloads/openbugs/OpenBUGS.zip

Now everything is reading to go. To launch WinBUGS:

wine ~/Downloads/WinBUGS14/WinBUGS14.exe

To launch OpenBUGS:

wine ~/Downloads/openbugs/winbugs.exe

Note the first line of code is for WinBUGS and the second line is for OpenBUGS.

Don’t forget to add the immortality key and the patch for WinBUGS.

h1

Bayesian Guide for Education

October 24, 2009

I have been working on a gentle, very applied introduction to Bayesian aimed at educational researchers. I’ve uploaded the guide. Feel free to download it and if you do please provide feedback.

Click here for the Bayesian guide

h1

R just got a little more love from Debian

July 30, 2009
Announcing cran2deb: 1700+ Debian packages from almost all of CRAN
------------------------------------------------------------------

Last Friday's presentation at UseR! 2009 was the first really public mention
of 'cran2deb'. It provides Debian packages of all of CRAN. It started as
Charles' project from last year's Google Summer of Code, was further extended
by us over the last year.

We currently build for the testing distribution and the i386 and amd64
architectures. This is now publically useable and we welcome wider testing
by Debian users. We hope to extend this to Ubuntu during the summer or fall.

I include what I blogged earlier (which you may have seen via Planet Debian
and/or Planet R).  This also contains the /etc/apt/sources.list entries:

   cran2deb: Would you like 1700+ new Debian / R packages ?

   As I mentioned in my [8]quick write-up of UseR 2009, one of my talks
   was about cran2deb: a system to turn (essentially) all [1]CRAN packages
   into directly apt-get-able binary packages.

   This is essentially a '2.0' version of earlier work with Steffen
   Moeller and David Vernazobres which we had presented in 2007. Then, the
   approach was top-down and monolithic which started to show its limits.
   This time, the idea was to borrow the successful bottom-up approach of
   my [2]CRANberries feed.

   The bulk of the work was done by Charles Blundell as part of his Google
   Summer of Code 2008 project which I had suggested and mentored. After
   that project had concluded, we both felt we should continue with it and
   bring it to 'production'. The CRAN hosts provided us with a (virtual
   Xen) machine to build on, and we are now ready to more publically
   announce the availability of the repositories for i386 and amd64:

      deb http://debian.cran.r-project.org/cran2deb/debian-i386 testing/

   and

      deb http://debian.cran.r-project.org/cran2deb/debian-amd64 testing/

   A few more details are provided in [3]our presentation slides. We look
   forward to hearing from folks using; the r-sig-debian list may be a
   good venue for this.

   References
   1. http://cran.r-project.org/
   2. http://dirk.eddelbuettel.com/cranberries
   3. http://dirk.eddelbuettel.com/papers/useR2009cran2deb.pdf

Please direct questions to r-sig-debian (which has subscriber-only posts) or
debian-science. 

Regards,  Dirk and Charles

--
Three out of two people have difficulties with fractions.

I would like to thank Dirk and Charles for their hard work!

Original message posted by Dirk:

https://stat.ethz.ch/pipermail/r-sig-debian/2009-July/000805.html

h1

Simple OpenBUGS MR example

June 17, 2009

Below is some OpenBUGS code that I created that models income as a function of age, race, and educational status. This is a simple MR example and can easily be replicated in R with the lm() function.

## OpenBUGS Model
model{
for( i in 1:6641) {
income[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 + beta1*educ[i] + beta2*age[i] + beta3*race[i]
}
## Priors
#beta0 ~ dflat()
#beta1 ~ dflat()
#beta2 ~ dflat()
#beta3 ~ dflat()
## Hypothesized Priors, deduced by "peaking"
beta0 ~ dnorm(3, 0.1)
beta1 ~ dnorm(0.3, 0.1)
beta2 ~ dnorm(0.01,0.5)
beta3 ~ dnorm(0, 0.1)
tau ~ dgamma(0.1,0.1)
sigma <- 1/sqrt(tau)
}
## Initial values
list(beta0 = 0, beta1 = 1, beta2 = 1, beta3 = 1, tau = 1)

You can get the data and a file that contains the model and the initial values here. This is example is entirely didactic, comes with no warranty, and is entirely intended to help get folks up and running in *BUGS.

The most interesting bits: With this large sample size (n ~ 6000), the priors seem to be washed away by the data and the high similarity between the results from the Bayesian and the frequentist analyses (run via lm in R).

h1

OpenBUGS and WinBUGS on Mac OS X

June 15, 2009

I haven’t found a really good how-to for getting OpenBUGS and WinBUGS up and running on Mac OS X, so I thought I’d create a little how-to. Setting up these programs is super easy but requires MacPorts and wine .

First, you need to have the Xcode developer tools installed. Go to Apple, create an account, and download them.

Second, install MacPorts (get it here for 10.5) by clicking the *.dmg file.

Third, open a terminal and type:

sudo port install wine

This will take about 20 minutes to compile wine and all its dependencies.

Fourth, download the WinBUGS zip file here and the OpenBUGS zip file here .

Fifth, unzip the WinBUGS file by double clicking on it and it will create a folder called WinBUGS14. Move this to Applications. For OpenBUGS, I had to create a folder called OpenBUGS, move the zip file there, open up a terminal, and type unzip OpenBUGS.zip (double clicking this file did not work thus I fired up a terminal). This folder I then moved to Applications.

Sixth, to start WinBUGS open a terminal and type the following:

cd /Applications/WinBUGS14
wine WinBUGS14.exe

For OpenBUGS, open a terminal and type the following (assuming you created a folder named OpenBUGS with all the file from the OpenBUGS.zip files in them):

cd /Applications/OpenBUGS
wine winbugs.exe

Everything seems to work fine. The only thing I’ve noticed is that the fonts are not as nice in wine as they are in Windows and in WinBUGS some of the documentation doesn’t work properly (e.g. clicking on links does nothing), whereas in OpenBUGS they work flawlessly.

Please leave comments if this works/doesn’t work for you.

h1

IES 2009 Poster Presentation

June 11, 2009

This past Monday I presented a poster on my research at the Institute of Education Sciences . It examined the influence of risk on reading and math achievement with a specific focus on treating risk as a dynamic and static covariate. Analysis was performed with the lme4 package in R and graphs were created in ggplot2 . Model selection was done with the BIC. The static models overwhelmingly beat the dynamic models. Click on the picture to download a PDF of the poster. Feedback is welcomed.

IES2009

Finally, this blog will now be updated at least weekly with more comprehensive notes on Bayesian analysis, R and WinBUGS code.

h1

Bayesian Hierarchical Modeling

May 15, 2009

I’ve now dug a bit deeper into Carlin & Louis (2008) compliments of the end of the semester and have finally gotten around to Bayesian hierarchical modeling. The extensions of Bayesian models to more levels is very easy to do in WinBUGS (or OpenBUGS). Really the only thing you need to specify are priors and your model. Specification of priors seems to be a relatively simple thing to do in most cases. Priors for most model encountered in social sciences seem often follow Gamma distributions. Adding random effects to a model just requires additional priors for these variables (I believe they’re referred to as hyperpriors). So in the model you’ll have priors for precision (the reciprocal of variance), the random effects (often Gamma or Normal), and the fixed effects (often Normal or you could specify a flat or uniform prior). It’s important to try and pick a conjugate prior. Beyond priors, you just need to specify your model and finally initial values and WinBUGS is ready to roll. Don’t forget you have to burn in before you simulate. The authors seem to always run MCMC 1000 times before actual simulations.

Later today or in the next few days I plan to run a hierarchical model with achievement as my outcome, students and schools as random effects, and risk status and gender as my independent variables. I will post the WinBUGS code and discuss a little of the model selection done with DIC. Also I will play around with prior specification and see how that affects my estimates and Bayesian CIs.  Note: All the data shared here is simulated and does not represent any projects that I am currently involved.

h1

Bayesian Longitudinal Modeling in R

May 1, 2009

A lot of my interests and the data I seem to get my hands on are longitudinal. I frequently use the lme4() package in R to do my analyses (with both discrete and continuous data). However, I recently learned about a package called MCMCglmm() written by Jarrod Hadfield, which allows a full Bayesian analysis on longitudinal discrete or continuous data. I strongly recommend this package and Jarrod is incredibly helpful when you encounter walls.

To use the package:
install.packages("MCMCglmm")
library(MCMCglmm)
my.model <- MCMCglmm(y ~ 1 + predictor + time + time*predictor, random= ~ id, family="poisson", data=mydata, verbose=FALSE, DIC=TRUE)

This code will model my dependent variable (y) by time, my predictor, and a time-predictor interaction. I have allowed only random intercepts in my model (random = ~ id), and specified the poisson family. The more crucial bit of code is the DIC=TRUE, if you want to use DIC for model comparison.

To plot the posterior distributions of my parameters:
plot(my.model$Sol)

To see the posterior mode of the parameters:
posterior.mode(my.model$Sol[,1] ## For Intercept
posterior.mode(my.model$Sol[,2] ## For Predictor (lets assume 2 level dummy)
posterior.mode(my.model$Sol[,3] ## For Time
posterior.mode(my.model$Sol[,4] ## For Time*Predictor

Credible intervals can be obtained by the following:
HPDinterval(my.model$Sol[,1], 0.95) ## For Intercept … and so on …

For the DIC:
my.model$DIC

For more information about this package read the vignette, e.g. specifying priors, examining random effects.
library(MCMCglmm)
vignette("Tutorial")