Hi:

On Fri, Aug 27, 2010 at 9:12 AM, Dave Schoeman <d.schoeman@ulster.ac.uk>wrote:

> I tried this in a slightly different context and got no response, so I
> thought I'd try and frame the analysis slightly differently. The design is
> quite simple, but the analysis is proving difficult. Say I have 20
> individual Subjects prone to mild anxiety attacks, half of which are
> assigned to a Control group and half to an Experimental group. I monitor
> these patients for 8 consecutive years, counting the number of anxiety
> attacks each year. After the initial 4 years, I provide members of the
> Control group a placebo and members of the Experimental group a supposed
> effective treatment. The hypothesis is that the supposed effective treatment
> has no effect in reality.
>
> As I see it, I have several fixed effects, all of which are factors: Group
> (Control or Experiment); BA (Before "treatment" or After "treatment"); and
> Year (nested within BA). There is a single random effect: Subject (nested
> within Group).
>
> Data can be simulated as follows (I have coded the attack rate Before
> slightly higher in the Experiment than the Control Group, with the attack
> rate increasing slightly after dosing in the Control group, but decreasing
> substantially in the Experiment group):
>        Subject <- c(rep(c(1:20), 8))
>        Group <- c(rep(c(rep("C", 10), rep("E", 10)), 8))
>        BA <- c(rep("Bef", 80), rep("After", 80))
>        Year <- c(rep(1990, 20), rep(1991, 20), rep(1992, 20), rep(1993,
> 20), rep(1994, 20), rep(1995, 20), rep(1996, 20), rep(1997, 20))
>        Attack <- c(rep(c(rpois(10, lambda = 10), rpois(10, lambda = 12)),
> 4), rep(c(rpois(10, lambda = 12), rpois(10, lambda = 3)), 4))
>        dat <- data.frame(Attack, Subject, Group, BA, Year)
>        dat$Subject <- as.factor(dat$Subject)
>        dat$Year <- as.factor(dat$Year)
> The model seems like it should be:
>        mod <- lmer(Attack ~ Group*BA + BA/Year + BA/Year:Group +
> Group/Subject:BA + (1|Group/Subject), data = dat, family = poisson)
>
> Group is assigned at the time of the intervention; everything before it is
baseline data. I could see comparisons between before/after trends at the
subject level that would depend on Group assignment, but the Group * BA
interaction is a bit of a problem, since no treatment is assigned in the
'Before' phase. The (1 | Group/Subject) term makes some sense (random
intercepts), but you need something to generate and compare random slopes
before and after treatment by Group and Subject. I'm not quite sure how that
would be done; I'm concerned with maintaining the quantitative values of
Year but estimating the slopes and intercepts by before/after treatment
assignment. There must be a standard way to do this with biomedical data,
but I don't have a lot of experience in that area. Hopefully, others will
have better ideas.

I don't see how Group is nested within BA (BA/Year:Group), especially when
your first term is Group * BA, which suggests they are crossed. Group is not
even assigned temporally until the After phase of BA, so it's clearly not
nested within the levels of BA. It is more accurate to say that it is at
best partially crossed, but even that is debatable - the groups are split in
the After phase. You can't have it both ways - either two terms are crossed
or they are nested...or in this case, perhaps somewhere in between. This is
the part you need to unravel properly. Conceptually, I can see why you think
Year should be nested within BA, but Year is a quantitative variable rather
than a factor, and your aim should be to predict random intercepts and
slopes among subjects, and compare mean differences between the control and
experimental groups, before and after the treatments were applied.

This type of problem must occur in biostatistics at least occasionally. How
is this traditionally handled when you have longitudinal data by subject,
where part of it is baseline (pre-treatment) and part of it is
post-treatment?  (BA * Year | Group/Subject) ??


BUT, this gives an error: "Error in mer_finalize(ans) : Downdated X'X is not
> positive definite, 1."
>

This partially explains why:

library(ggplot2)
g <- ggplot(dat, aes(x = Year, y = Attack, groups = Subject))
g + geom_line()
g + geom_line() + facet_grid(Group ~ BA, scales = 'free')

You have no within-subject variablility within before/after phases and
several subjects overlap in at least one of the phases.

To avoid this error, I have to drop all of the interaction terms involving
> nesting. HOWEVER, if I code Year as a continuous covariate, I need only drop
> the term Group/Subject:BA. This suggests that interaction terms comprising
> only factors are not easily fit using this unreplicated design...
>

I believe you have hit the nail on the head....rather than unreplicated, I
would say incompletely crossed [BA * Group], because there is replication in
the BA and Group factors.

>
> My questions are these:
> 1 - Does my initial model look correct (and do I need to specify the
> nesting - other threads suggest that this might not be necessary if nested
> variables have unique codes)?
> 2 - Why is it that the  factor-by-factor interactions won't work (I'm an
> lmer() novice, so it might be very obvious to some...)?
> 3 - Would I be more or less correct to simply ignore the Group/Subject:BA
> interaction and proceed with model development from there?
>
> Any comments appreciated.
>
> - Dave
>
>
I don't think I've helped resolve the form of model you need, but it seems
that certain aspects of your design have to be massaged. Hopefully others
can weigh in on the problem...

Regards,
Dennis

>
>
> _________
> Dr Dave Schoeman
> Lecturer in Marine Science
> School of Environmental Sciences
> University of Ulster
> Cromore Road
> Coleraine
> BT52 1SA
> Northern Ireland
>
> Phone: +44 (0)28 701 24076
> Mobile: +44 (0)75 49 526 743
> Fax: +44 (0)28 701 24491
>
> Treat data with care; under duress they will tell you anything you want to
> hear...
>
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]


