[R-sig-ME] nAGQ = 0

[Rolf Turner]

On 04/09/17 10:51, Poe, John wrote:

<SNIP>

> You can try using the BRMS package if you aren't comfortable switching 
> to something totally unfamiliar. It's a wrapper for Stan designed to use 
> lme4 syntax and a lot of good default settings. It's pretty easy to use 
> if you know lme4 syntax and can read up on mcmc diagnostics.

<SNIP>

I remember a quote from some Brit comedy series (it starred the bloke 
who played Terry on "Minder"):  "My mother always said that you should 
try anything once, except for incest and Morris dancing.  She was right 
about the Morris dancing."

On that basis I decided to try using brms.  I installed the package from 
CRAN with no real problems, although there were some slightly worrying 
and highly technical warnings from "rstan" --- I think:

> warning: non-static data member initializers only available with -std=c++11 or -std=gnu++11
>        bool has_var_ = false;

Then I read the vignette brms_overview a bit, and plunged in with trying 
to fit a model.  Needless to say, the attempt didn't get much past 
square zero.  I tried it again with my artificial simulated data that I 
talked about in the post that started off this train of craziness (it 
elicited the suggestion from Tony Ives that I try using nAGQ=0).  The 
result was the same --- square zero + epsilon:

> Compiling the C++ model
> Start sampling
> 
> SAMPLING FOR MODEL 'binomial(cloglog) brms-model' NOW (CHAIN 1).
> Rejecting initial value:
>   Log probability evaluates to log(0), i.e. negative infinity.
>   Stan can't start sampling from this initial value.
     .
     .
     .
> Rejecting initial value:
>   Log probability evaluates to log(0), i.e. negative infinity.
>   Stan can't start sampling from this initial value.
> Rejecting initial value:
>   Log probability evaluates to log(0), i.e. negative infinity.
>   Stan can't start sampling from this initial value.
> 
> Initialization between (-2, 2) failed after 100 attempts. 
>  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
> [1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
> error occurred during calling the sampler; sampling not done

Since I'm flying completely blind here (no idea WTF I'm doing) I have 
come to a shuddering halt.

I have attached my function "artSim.R" to generate the artificial data,
and a script to source to effect the call to brm() that I used.

If some kind mixed models guru could take a look and point out to me 
just what bit of egregious stupidity I am committing, I'd be ever so 
humbly grateful.  I don't know from Bayesian stuff (priors, and like 
that) at all, so it's likely to be something pretty stupid and pretty 
simple in the first instance.

cheers,

Rolf Turner

-- 
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276
-------------- next part --------------
artSim <- function(){
#
# Function to simulate "artificial" data which is at least superficially
# similar to some real data.
#
require(MASS)

link    <- "cloglog"
B       <- binomial(link=link)
linkfun <- B$linkfun
linkinv <- B$linkinv

# Construct (artificial) treatment factor, covariate, and
# (random) replicate factor.
x    <- seq(0,28,by=2)
Trt  <- LETTERS[1:24]
Rep  <- 1:3 # Three reps per treatment.
Xdat <- expand.grid(x=x,Trt=Trt,Rep=Rep)
uRep <- with(Xdat,factor(paste0(Rep,Trt)))
Xdat$Rep <- with(Xdat,factor(as.numeric(uRep)))

beta0 <- seq(-3,0.45,by=0.15)
beta1 <- rep(seq(0.05,0.3,by=0.05),4)
names(beta0) <- Trt
names(beta1) <- Trt
Sigma <- matrix(c(0.06,-0.001,-0.001,0.0001),nrow=2)

lb0   <- beta0[match(Xdat$Trt,names(beta0))]
lb1   <- beta1[match(Xdat$Trt,names(beta1))]
nrep  <- 72
imat  <- match(Xdat$Rep,1:nrep)
Z     <- mvrnorm(nrep,c(0,0),Sigma)[imat,]
linpr <- lb0 + Z[,1] + (lb1 + Z[,2])*Xdat$x
p     <- linkinv(linpr)
nsize <- 25
Dead  <- rbinom(nrow(Xdat),nsize,p)
Alive <- nsize - Dead
x0    <- (linkfun(0.99) - beta0)/beta1
Xdat$Dead  <- Dead
Xdat$Alive <- Alive
attr(Xdat,"trueLD99") <- x0
return(Xdat)
}
-------------- next part --------------
#
# Script scr.brmsDemo
#

set.seed(42)
X <- artSim()
library(brms)
fit <- brm(Dead|trials(Dead+Alive) ~ (0+Trt)/x + (x | Rep),
           family=binomial(link="cloglog"),data=X,
           prior=set_prior("normal(0,10)"))