nAGQ = 0
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)"))