The rats data are described in the RatWeightData notebook where they are converted from a matrix to a saved R data frame. We could do the initial data manipulation in R through the RCall package for Julia but we chose to use packages from the "Hadleyverse" that cannot easily be installed on (https://juliabox.org)
In [1]:
using DataFrames, MixedModels, Mamba, RCall, Stan
In [2]:
versioninfo(true)
We retrieve the data as an R object, and copy it under the same name into Julia.
In [4]:
R"""
rats <- readRDS("./rats.rds")
""";
@rget rats
Out[4]:
At this point we do something radical and plot the data. I have never seen a data plot in any of the MCMC exampes that use these data. The panels are ordered (bottom to top, left to right) according to increasing average weight of the rat. The aspect ratio is chosen so that a typical slope of the within-rat least squares line is approximately 45 degrees on the plotting surface.
In [5]:
R"""
library(lattice)
print(xyplot(y ~ day | reorder(id, y), rats, type = c('p','g','r'),
aspect = 'xy', ylab = "Weight (g.)"))
""";
There is an overall linear trend in the weight with respect to time but there is also noticeable downward curvature for many of the rats. Nevertheless we will start with a model with linear model with vector-valued random effects for slope and intercept by rat.
The simplest way to write the model is weight ~ 1 + day + (1 + day|id) which allows for correlated random effects for slope and intercept for each rat.
In [6]:
R"""
suppressPackageStartupMessages(library(lme4))
m1 <- lmer(y ~ 1 + day + (1 + day | id), rats, REML = FALSE)
summary(m1)
"""
Out[6]:
Because day is days past birth, the estimate of a typical birth weight is 106.6 g. and the typical weight gain per day after birth is 6.2 g./day. As is common in linear regression models where x = 0 is to the left of the observed data, there is a negative correlation, -0.343, between these estimates. The standard deviation of the random effects for the intercept (i.e. the birth weight) is 10.5 g. and the standard deviation of the random effects for the slope is 0.50 g./day. There is a slight negative within-rat correlation, -0.15, of these random effects. We can check the conditional means of these random effects
In [7]:
rcall(:ranef, :m1)
Out[7]:
or, a better choice, plot these conditional modes.
In [8]:
R"""
re1 <- ranef(m1, condVar = TRUE)
print(dotplot(re1, scales = list(x = list(relation = 'free'))))[[1]]
""";
We can check for non-negligible correlation of the random effects by fitting a model with uncorrelated random effects and comparing the fits.
In [9]:
R"""
m2 <- lmer(y ~ 1 + day + (1|id) + (0 + day|id), rats, REML = FALSE)
options(show.signif.stars = FALSE)
anova(m2, m1)
"""
Out[9]:
In [10]:
rcall(:summary, :m2)
Out[10]:
In [11]:
m1 = fit!(lmm(y ~ 1 + day + (1 + day | id), rats))
Out[11]:
The two mixed-models packages agree on the fit. The amount of time required for the lmm fit is
In [12]:
@time fit!(lmm(y ~ 1 + day + (1 + day | id), rats));
In [13]:
m2 = fit!(lmm(y ~ 1 + day + (1 | id) + (0 + day | id), rats))
Out[13]:
In [14]:
MixedModels.lrt(m2, m1) # print format is still a bit primitive
Out[14]:
The independent random effects model, m2, has been a standard example for MCMC methods since the original BUGS system. The Stan examples respository, https://github.com/stan-dev/example-models, contains links to this example in OpenBUGS and in Stan, under bugs-examples/vol1/rats/.
The Stan package for Julia is similar to rstan in that the interactive language is used to marshall the data which is then passed to a stan program. The initial version of the model uses the data in the form of a matrix of response values, y, and a separate vector of times called x. Because Stan is a statically-typed language (it is transformed to C++ code), the sizes of all arrays must be given explicitly. The data should be available as a Dict{ASCIIString,Any} type.
In [14]:
const ratsdict = Dict(
"N" => 30,
"T" => 5,
"x" => [8.,15,22,29,36],
"xbar" => 22.,
"y" => reshape(convert(Vector{Float64}, rats[:y]), (30, 5)))
Out[14]:
In [15]:
ratstan = Stanmodel(name = "rats", model = """
# http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/Vol1.pdf
# Page 3: Rats
data {
int<lower=0> N;
int<lower=0> T;
real x[T];
real y[N,T];
real xbar;
}
parameters {
real alpha[N];
real beta[N];
real mu_alpha;
real mu_beta; // beta.c in original bugs model
real<lower=0> sigmasq_y;
real<lower=0> sigmasq_alpha;
real<lower=0> sigmasq_beta;
}
transformed parameters {
real<lower=0> sigma_y; // sigma in original bugs model
real<lower=0> sigma_alpha;
real<lower=0> sigma_beta;
sigma_y <- sqrt(sigmasq_y);
sigma_alpha <- sqrt(sigmasq_alpha);
sigma_beta <- sqrt(sigmasq_beta);
}
model {
mu_alpha ~ normal(0, 100);
mu_beta ~ normal(0, 100);
sigmasq_y ~ inv_gamma(0.001, 0.001);
sigmasq_alpha ~ inv_gamma(0.001, 0.001);
sigmasq_beta ~ inv_gamma(0.001, 0.001);
alpha ~ normal(mu_alpha, sigma_alpha); // vectorized
beta ~ normal(mu_beta, sigma_beta); // vectorizedsim1 = stan(ratstan);
for (n in 1:N)
for (t in 1:T)
y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), sigma_y);
}
generated quantities {
real alpha0;
alpha0 <- mu_alpha - xbar * mu_beta;
}
""",
data = [ratsdict],
monitors = ["mu_alpha", "mu_beta", "sigma_y", "sigma_alpha", "sigma_beta"])
Out[15]:
In [16]:
display(ratsdict)
In [17]:
sim1 = stan(ratstan, ratsdict)
Out[17]:
In [18]:
## Data
rats1 = Dict{Symbol, Any}(
:y => convert(Vector{Float64}, rats[:y]),
:rat => convert(Vector{Int}, rats[:id]),
:X => convert(Vector, rats[:day])
)
rats1[:xbar] = mean(rats1[:X])
rats1[:Xm] = rats1[:X] - rats1[:xbar]
display(rats1)
In [19]:
## Model Specification
model = Model(
y = Stochastic(1,
(alpha, beta, rat, Xm, s2_c) ->
begin
mu = alpha[rat] + beta[rat] .* Xm
MvNormal(mu, sqrt(s2_c))
end,
false
),
alpha = Stochastic(1,
(mu_alpha, s2_alpha) -> Normal(mu_alpha, sqrt(s2_alpha)),
false
),
alpha0 = Logical(
(mu_alpha, xbar, mu_beta) -> mu_alpha - xbar * mu_beta
),
mu_alpha = Stochastic(
() -> Normal(0.0, 1000),
false
),
s2_alpha = Stochastic(
() -> InverseGamma(0.001, 0.001),
false
),
beta = Stochastic(1,
(mu_beta, s2_beta) -> Normal(mu_beta, sqrt(s2_beta)),
false
),
mu_beta = Stochastic(
() -> Normal(0.0, 1000)
),
s2_beta = Stochastic(
() -> InverseGamma(0.001, 0.001),
false
),
s2_c = Stochastic(
() -> InverseGamma(0.001, 0.001)
)
);
In [20]:
inits = [
Dict(:y => rats1[:y], :alpha => fill(250, 30), :beta => fill(6, 30),
:mu_alpha => 150, :mu_beta => 10, :s2_c => 1, :s2_alpha => 1,
:s2_beta => 1),
Dict(:y => rats1[:y], :alpha => fill(20, 30), :beta => fill(0.6, 30),
:mu_alpha => 15, :mu_beta => 1, :s2_c => 10, :s2_alpha => 10,
:s2_beta => 10)
];
In [21]:
scheme = [Slice(:s2_c, 10.0),
AMWG(:alpha, 100.0),
Slice([:mu_alpha, :s2_alpha], [100.0, 10.0], Univariate),
AMWG(:beta, 1.0),
Slice([:mu_beta, :s2_beta], 1.0, Univariate)]
setsamplers!(model, scheme)
Out[21]:
In [22]:
sim = mcmc(model, rats1, inits, 10000, burnin=2500, thin=2, chains=2);
describe(sim)