[R-sig-ME] nAGQ = 0
Rolf Turner
r.turner at auckland.ac.nz
Wed Sep 6 23:48:08 CEST 2017
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)"))
More information about the R-sig-mixed-models
mailing list