
Thanks, Alain, for the advice and confirmation that I do need to use BUGS. 

For reference, I've had success using JAGS (R2jags) to run the WinBUGS code in your 2012 book. All the code I've tried runs fine as quoted using jags{R2jags} (the only exception so far is calculating Pearson residuals: https://stat.ethz.ch/pipermail/r-sig-mixed-models/2012q3/019087.html)

For my current model, I've managed to get it running in JAGS without the temporal correlation. To add the AR correlation I'm using your example in Section 9.8.2 of your 2012 book. The book code runs fine in JAGS, but as soon as I try to convert to a binary response variable (converted Y to 0/1) by replacing 'log(mu[i])' with 'logit(mu[i])' and Y[i] ~ dpois(mu[i]) with Y[i] ~ dbin(mu[i],1), I get errors:  

Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains,  : 
  LOGIC ERROR:
Invalid mask in LogicalNode::isClosed

I know this is a jags error, but I suspect the problem is due to my model specification. Apologies if I'm making a beginner's mistake, but is there anything else I need to do to switch from a Poisson to binary response variable? The example works with Poisson response and AR, and with a binary response and no AR correlation. [I haven't got as far as the multinomial (ordered) response yet, but any tips would be greatly appreciated as examples seem to be a bit scarce!]

Thanks very much,

Andrew

####### Zuur et al. 2012 Ch9.8.2 example, with binary response:

require(R2jags)

# Read data
SW <- read.table("http://www.highstat.com/Book4/Spermwhales.txt", header=T, sep="\t")
MyVar <- c("Year", "Nsea", "AnnualTTAnom", "WinterNAO")
SW1 <- SW[SW$Year >=1563 & SW$Year <= 2000, MyVar]
SW1<-SW1[1:dim(SW1)[1]-1,MyVar]  # take off last row to avoid NA
SW1$TempAn <- as.numeric(scale(SW1$AnnualTTAnom))
SW1$NAO <- as.numeric(scale(SW1$WinterNAO))
N <- nrow(SW1)
num.knots <- 9
Prob <- seq(0, 1, length= (num.knots+2))
Prob1 <- Prob[-c(1, (num.knots+2))]
X <- cbind(rep(1,N), SW1$TempAn)
Knots <- quantile(unique(SW1$TempAn), probs=Prob1, na.rm=F)
GetZ <- function(x, Knots){
	Z_K <- (abs(outer(x, Knots, "-")))^3
	OMEGA_all <- (abs(outer(Knots, Knots, "-")))^3
	svd.OMEGA_all <- svd(OMEGA_all)
	sqrt.OMEGA_all <- t(svd.OMEGA_all$v %*% (t(svd.OMEGA_all$u) * sqrt(svd.OMEGA_all$d)))
	Z <- t(solve(sqrt.OMEGA_all, t(Z_K)))
	return(Z)
}
Z <- GetZ(SW1$TempAn, Knots)

win.data <- list(Y = SW1$Nsea, NumKnots=9, X=X, Z=Z, N=nrow(SW1))
# Make response variable binomial
win.data2<-win.data
win.data2$Y[win.data2$Y > 0] <-1
sink("zuur2012.9-8-2.binar.bugs.txt")
cat("
    model{
      # Priors for regression parameters
      for (i in 1:2) {beta[i] ~ dnorm(0, 0.001)}
      for (i in 1:NumKnots) {b[i] ~ dnorm(0, taub)}
      taub ~ dgamma(0.001, 0.001)
      taueps ~ dgamma(0.001, 0.001)
      rho ~ dunif(-0.99, 0.99)
      
      # First row (i=1)
      # Calculate X*beta
      eta[1] <- beta[1] * X[1,1] + Smooth1[1]
      Smooth1[1] <- beta[2] * X[1,2] + b[1]*Z[1,1] + b[2]*Z[1,2] + b[3]*Z[1,3] + b[4]*Z[1,4] + b[5]*Z[1,5] * b[6]*Z[1,6] * b[7]*Z[1,7] * b[8]*Z[1,8] * b[9]*Z[1,9]
      # Calculate X*beta + eps 
      eta1[1] ~ dnorm(eta[1], taueps)
      logit(mu[1]) <- eta1[1]
#      Y[1] ~ dbern(mu[1])
      Y[1] ~ dbin(mu[1],1)
    
      for(i in 2:N){
      eta[i] <- beta[1] * X[i,1] + Smooth1[i]
      Smooth1[i] <- beta[2] * X[i,2] + b[1]*Z[i,1] + b[2]*Z[i,2] + b[3]*Z[i,3] + b[4]*Z[i,4] + b[5]*Z[i,5] * b[6]*Z[i,6] * b[7]*Z[i,7] * b[8]*Z[i,8] * b[9]*Z[i,9]
        # Calculate eta + rho*(logit(u[i=1]) - eta[i-1]) + eps
        temp[i] <- eta[i] + rho * (logit(mu[i-1]) - eta[i-1])
        eta1[i] ~ dnorm(temp[i], taueps)
        logit(mu[i]) <- max(-10, min(10, eta1[i]))
#        Y[i] ~ dbern(mu[i])
        Y[i] ~ dbin(mu[i],1)

      }
    }
    ", fill=TRUE)
sink()

inits <- function(){
	list( beta = rnorm(2, 0, 0.001), b = rnorm(9, 0, 0.001), taub=runif(1, 0.1, 0.5), taueps=runif(1, 0.1, 5), rho=runif(1, -0.99, 0.99))
}
params <- c("beta", "b", "Smooth1", "rho")

# Short run to test
mjags.binar <- jags(data=win.data2, inits=inits, parameters.to.save=params, model.file="zuur2012.9-8-2.binar.bugs.txt", n.chains=3, n.burnin=500, n.iter=3000, n.thin=1)

Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains,  : 
  LOGIC ERROR:
Invalid mask in LogicalNode::isClosed






[R-sig-ME] Is there an R function for GLMM with binary, response, nested random factors, and temporal correlation?

Highland Statistics Ltd highstat at highstat.com 
Wed Aug 29 12:26:46 CEST 2012

Previous message: [R-sig-ME] Fwd: MCMCglmm poisson / not poisson
Next message: [R-sig-ME] Interpretation of effects
Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]


>
> ------------------------------
>
> Message: 2
> Date: Tue, 28 Aug 2012 18:47:25 +1200
> From: Andrew Digby <andrewdigby at mac.com>
> To: Douglas Bates <bates at stat.wisc.edu>,
> 	"r-sig-mixed-models at r-project.org" <r-sig-mixed-models at r-project.org>
> Subject: Re: [R-sig-ME] Is there an R function for GLMM with binary
> 	response, nested random factors, and temporal correlation?
> Message-ID: <14D7DCCA-01B9-433E-9727-672B50E559DB at mac.com>
> Content-Type: text/plain; CHARSET=US-ASCII
>
>
> Thanks again, and apologies for the delay. I'd prefer not to incorporate the temporal correlation structure through a random intercept/slope as I want to try auto-regressive correlation rather than compound symmetry.
>
> I've tried MCMCglmm, but it doesn't seem to allow AR correlation either. I don't have time to get into running WinBugs on OS X, or

After 5 days fulltime pottering around last year, I gave up trying to 
get WinBUSG to run on OS X (actually...I did manage to run OpenBUGS on 
OS X). Now I use TeamViewer (free)...which allows me to remotely access 
a Windows machine from a Mac, and run for example WinBUGS. It takes 
about 1 minute to set up. And WinBUGS runs through the night.

All you need is remote access to an Windows computer...which can't be 
too difficult to arrange.



> learning ADMB. Any further suggestions on how to do mixed logistic regression with temporal correlation would be very welcome!
I'm afraid you will have to go the WinBUGS/OpenBUGS route. Example code 
for what you want is in our 2012 book.

Kind regards,

Alain


> Thanks,
>
> Andrew
>
> On 10/08/2012, at 04:21 , Douglas Bates <bates at stat.wisc.edu> wrote:





>> On Wed, Aug 8, 2012 at 7:49 PM, Andrew Digby <andrewdigby at mac.com> wrote:
>>> Thank you for your advice. Can I interpret it  to mean that it's not correct to specify a correlation structure (beyond that induced through random effects) in generalised models? I had thought that although prone to some issues, this was possible; see e.g. Ben Bolker's and Alain Zuur's reponses here: https://stat.ethz.ch/pipermail/r-sig-mixed-models/2012q1/014937.html. In my case, I think I do have 'small-scale correlation', with repeated measures of a subject taken in rapid succession on a particular day. My residuals don't show gross heterogeneity, so I may not need to alter the variance structure to allow for that.
>> Perhaps Ben or Alain can provide more information on models that
>> incorporate correlation structure in binary responses.  There may be
>> definitions of these using copulas or something like that but I don't
>> know enough about that area to be able to provide advice.
>>
>>> I'd be grateful for general advice on whether a GLM with temporal correlation structure is an adequate method in this case. Or should I tackle MCMC instead? (I have Zuur et al's excellent books.) However, this is a small part of my overall analysis, so I'd prefer to keep things as simple as possible!
>> One way to achieve a temporal correlation structure is to have a
>> random effect for slope with respect to time for each subject.
>>
>>> On 9/08/2012, at 03:19 , Douglas Bates <bates at stat.wisc.edu> wrote:
>>>
>>>> On Tue, Aug 7, 2012 at 5:21 PM, Andrew Digby <andrewdigby at mac.com> wrote:
>>>>> Despite lots of investigation, I haven't found any R packages might be suitable for the following problem. I'd be very grateful for suggestions.
>>>>>
>>>>> I have three-way nested data, with a series of measures (obs) taken in quick succession (equal time spacing) from each subject on different days. The measures taken on the same day are temporally correlated, so I'd like to use an AR1 correlation structure for those, but treat subjects and days as nested random factors (random intercept) since there is little temporal correlation between days. The response is binary.
>>>>>
>>>>> So I need a GLMM with a correlation structure. I've tried using GEE, but the R packages can't cope with multilevel nested data. The only R function I've found that can do this is glmmPQL.
>>>> Before you look for an R function, you should first check whether
>>>> there is indeed a statistical model with the properties that you
>>>> mention.  In the standard definition of a generalized linear mixed
>>>> model, and the only one that makes sense to me, the conditional
>>>> distribution of the response given the random effects has independent
>>>> components.
>>>>
>>>> Aspects of linear mixed models that depend on being able to model the
>>>> variance-covariance of the response separately from the mean don't
>>>> carry over to generalized linear mixed models.  One of the fundamental
>>>> properties of GLMs and GLMMs is that the variance does depend on the
>>>> mean.
>>>>
>>>>> m <- glmmPQL(y ~ f1 * f2 * f3 + (1|subj/day), correlation=corAR1(form =~obsno|subj/day))
>>>>>
>>>>> f1 - f3 are fixed factors
>>>>>
>>>>> However, PQL estimation is not recommended for binary response data. With no AIC and unreliable p values, model selection seems impossible! So my question is:
>>>>>
>>>>> 1) are there any other functions which are suitable for a GLMM with multilevel nested random effects and a AR1 correlation structure? Or is MCMC the only option?
>>>>> 2) to make things more complicated, I'd also like to include a varFunc variance structure to cope with heterogeneity. Is this possible in ML methods in R? I'd also like to extend to a multinomial response at a later stage.
>>>>>
>>>>> GEE seems the best bet, but I come unstuck with the three-way nested factors.
>>>>>
>>>>> Thanks for your help.
>>>>>
>>>>> Note: I originally posted this on R-help, but it was suggested that this list might be more appropriate.
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-mixed-models at r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>




********************************





-- 

Dr. Alain F. Zuur
First author of:

1. Analysing Ecological Data (2007).
Zuur, AF, Ieno, EN and Smith, GM. Springer. 680 p.
URL: www.springer.com/0-387-45967-7


2. Mixed effects models and extensions in ecology with R. (2009).
Zuur, AF, Ieno, EN, Walker, N, Saveliev, AA, and Smith, GM. Springer.
http://www.springer.com/life+sci/ecology/book/978-0-387-87457-9


3. A Beginner's Guide to R (2009).
Zuur, AF, Ieno, EN, Meesters, EHWG. Springer
http://www.springer.com/statistics/computational/book/978-0-387-93836-3


4. Zero Inflated Models and Generalized Linear Mixed Models with R. (2012) Zuur, Saveliev, Ieno.
http://www.highstat.com/book4.htm

Other books: http://www.highstat.com/books.htm


Statistical consultancy, courses, data analysis and software
Highland Statistics Ltd.
6 Laverock road
UK - AB41 6FN Newburgh
Tel: 0044 1358 788177
Email: highstat at highstat.com
URL: www.highstat.com
URL: www.brodgar.com


Previous message: [R-sig-ME] Fwd: MCMCglmm poisson / not poisson
Next message: [R-sig-ME] Interpretation of effects
Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
More information about the R-sig-mixed-models mailing list

	[[alternative HTML version deleted]]

