[R-sig-ME] Two Questions re: Piecewise Mixed Effects Models

[Joshua Wiley]

Hi Eric,

The first part of this code is reproducible and I suggest you run it
to see the graph and the output as I hope it will help you understand.

## simulate some data
set.seed(10)
y <- rnorm(1000, mean = (time <- sample(-2:2, 1000, TRUE)) * ((time >
0) + .5), sd = 1)

# plot time and the outcome y, you can see the 'break' sort of
plot(jitter(time), y)

# fit a non broek and broken time model
# also get the predicted values for times -2, -1, 0, 1, and 2
yhat1 <- predict(m1 <- lm(y ~ time), data.frame(time = -2:2))
yhat2 <- predict(m2 <- lm(y ~ time + time:I(time > 0)), data.frame(time = -2:2))

# add lines of predicted values from the different models
lines(-2:2, yhat1, col = "black", lwd = 3)
lines(-2:2, yhat2, col = "blue", lwd = 3)

# summaries of the models so you can see the coefficients
summary(m1)
summary(m2)

# what you have is an intercept and slope for time <= 0
# these are called intercept and time
# then for time > 0, you get an adjustment to the time slope (the
interaction term)
# the final design matrix would look something like:
cbind(Int = 1, Time = -2:2, Time2Adj = (-2:2) * (-2:2 > 0))


This part of the code is not reproducible, but I think this is more
like the model you probably want in glmer()

model <- glmer(rate ~ policy + time.before + time.after + Z +
  (1+ time.before + time.after | Site) + offset(log(people)),
  data=data, family="poisson")

# this gives you fixed effects for policy and Z
# and random effects for the intercept, and two time slopes
# that are allowed to vary across levels in Site
# all random effects have freely estimates covariances, so that is a 3
x 3 matrix
# (intercept, time.before, time.after) with all elements estimated

This is roughly (exactly?) equivalenet to the SAS model using proc glimmix:

proc glimmix data=data method=LAPLACE;
  class Site;
  model rate = policy time.before time.after Z / dist=poisson offset=people;
  random intercept time.before time.after / subject=Site type = chol;
/* you may be more familiar with type = un for unstructured */
/*  random intercept time.before time.after / subject=Site type = un; */
run;

Hope this helps,

Josh


On Sun, Sep 30, 2012 at 2:20 AM, Lofgren, Eric <elofgren at email.unc.edu> wrote:
> I'm relatively new to mixed effects models, and have two fairly basic questions (I think) about a project I'm working on.
>
> Essentially, I'm working on a poisson regression problem where we're trying to estimate the impact of a change in a policy on the rate of an event X. There are N sites where this change has taken place, and we also happen to know another characteristic of each site Z. It's been suggested that I look a a piecewise mixed effects model (which I've also seen referred to as an interrupted time series or "broken stick" model).
>
> I've got two questions about this type of model:
>
> 1. As I understand it, you should have two time variables, a "Time Before" and a "Time After" variable, each of which equals 0 right at the point of the policy change. My question is what these variables should look like, specifically what they should look like after they reach zero. For example, with a simple series of 5 points and the policy change in the middle, should it look like this:
>
> Time Before: 2, 1, 0, NA, NA
> Time After: NA, NA, 0, 1, 2
>
> or this?
>
> Time Before: 2, 1, 0, -1, -2
> TIme After: -2, -1, 0, 1, 2
>
> Or something else entirely?
>
> 2. In terms of using something like glmer() to fit the model itself, I'm a little confused about the syntax (I normally use SAS, but am trying to get better at using R). Following some advice I've seen, I should be using something like this:
>
> model <- glmer( rate ~ policy + time.before + time.after + (1+ time.before + time.after + Z |Site) + offset(log(people)), data=data, family="poisson")
>
> First, is this correct? Second, what's accomplished by having the time variables in both the fixed and random effects sections - what does that produce? A fixed slope and random intercept? Something else?
>
> Thanks,
>
> Eric
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/