Author Archives: nrc697sa-finnj

Simulation of More Complex Models

Hi class,

The book by Marc Kery “Introduction to WinBUGS for Ecologists” makes extensive use of simulated datasets. It starts with a mean model, and then t-tests, and goes up through GLMMs including ZIPs and mixture models.   See the book web site for a couple of sample chapters.

The code from his book is in this file R_WB_Code.txt .  There are no data files because he simulates everything (I think he does use real data in some of the exercises).  I also downloaded an errata file and an updated appendix with tips for using WinBUGS.  WinBUGS and OpenBUGS are very similar, so most of the advice should apply to either program.

I’ll go over one of the examples for simulating a GLMM.  We’ll even throw in some autocorrelation, if I can figure out how to do it.

Jack

Example MCMC Using OpenBUGS

Hi class,

Here is an example script following Zuur’s Chapter 1 from his new book ( Chap1_jtf.r ).  I added a few things, including using the coda package to plot some of the BUGS output.  The data on feeding of oyster catchers is in the file OystercatcherData.txt (from Zuur’s website http://highstat.com .  Also, here is an  Appendix from “Bayesian Modeling  Using WinBUGS,” by Ioannis Ntzoufras, on how to use the R package “coda”.

Enjoy your weekend!

Jack

Simulating A Fake Dataset

Hi class,

Here is the script I played with today: Chapter11_jtf.r

I used the M1 model to simulate a couple of fake datasets.

Gelman and Hill spend a fair amount of space in their book discussing simulation, and they have some software to do it.

As random structures get more complicated, simulation becomes harder, but I think we can steal it from them.  More on this later.

Jack

Logistic Regression on Tree Failure – Weird Results

Hi class,

I’ve been working on Brian Kane’s tree failure data from last October’s snowstorm, and got some weird results.  I set up a logistic regression with Failure as a response (tree broke = 1, not=0), and something funny is happening.  When I run a full model with a random slope and intercept, it doesn’t run on my portable, but it does run on my desktop!  Attached are the data file for Acer rubrum: 190 trees in several locations around campus.  The locations are the random effect.  DBH, presence of a defect and leaves on or off are the three fixed effects.  Here is my data file (dfAcRu.csv) and my R script (AcerRubrum.r).  If you have a chance, could you try running the script on your computer and see if it works?

Thanks.

Jack

Bates’s Book on Mixed Effects Models

Hi class,

The author of the lme4 package has an online book that is worth taking a look at.  The latest and greatest version resides in the R-forge web site, and you should load lme4 from there!  Some code and the chapters he has put up on the R-forge site are in the zip document Bates_MixedEffectsModels.zip .  We’ll talk about the book and some new commands on Tuesday.

Jack

 

Bolker’s Reanalysis of Owl Data

Hi class,

Here is an interesting example from Ben Bolker using the owls dataset.   He shows how to use an observation level group effect (to effectively change to a Poisson-lognormal distribution), how to predict the expected values of a lme4 model, and an interesting way to plot expected vs. observed with confidence intervals.  I might have to learn ggplot2! Enjoy!

Here are the pdf and the script I did in class:

Owls.pdf

BolkersOwls.r

Jack

Chapter 13 GLMM and GAMM

Hi class,

Here is my modified script for Chapter 13 – Chapter13_jtf.r

Can’t quite get to everything we might want to do.  Here is a paper that talks about GLMMs: Bolker et al. 2009 .  It goes into some detail about how each pacakge (or program) solves the problem.  Apparently, there is real difficulty in solving GLMMs and GAMMs.  The difficulty seems to be that the solution to the fixed effects have to be integrated over the random effects.  This is not exact, but there are various approximations.  PQL (penalized quasi-likelihood) is worst, LaPlace Approximation is better, etc.

Also, look at this link (http://glmm.wikidot.com/pkg-comparison#GLMMs), which lists the various commands and what you can do with them.  For example, it looks like glmmPQL is the only command that can do correlation structures.

Have fun!

Jack

Chapter 12 – Revised Script

Hi class,

Here is the final script for Chapter 12 (Chapter12_jtf_v2.r).  I made a few changes in class and after class, I figured out why I was confused about the QIC and yags.  It looks like yags (using QIC) and geeglm (using anova) come to different conclusions about whether to drop the interaction. The QIC says don’t, and the anova says do.  It’s hard to know if the problem is the software (yags vs geeglm) or the metric (likelihood ratio vs. QIC).

I think the bottom line on GEE is that it is very useful to try, especially if you can’t get just the right distribution to work.  It can mimic mixed effects (using intraclass correlation), can account for various other correlation structures (like AR1) and may be able to give you most of what you want.  I’d stick to the anova approach, unless you want to look at the references in detail.

Jack