Hi John,
It sounds like what you need is a mixed effects model. This will allow you to
model effects which apply to the whole population (fixed effects) and those
which are specific to an individual (random effects).The theory and practice of
mixed effects models is too large to discuss fully in an email. But I can point
you in the right direction.
First of all I suggest you read "Mixed-Effects models in S and S-PLUS" by
Pinheiro and Bates. This is an excellent text on the subject with loads of
worked examples. Also John Fox has produced a very nice introduction to mixed
effects models, again with worked examples in R: see
http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-mixed-models.pdf
The library you require in R is called library(nlme).
look up the Pixel example (?Pixel) it shows you how to model non-linear effects
using a quadratic term.
Download the nlme library from CRAN and paste the following in to R:
library(nlme)
fm1 <- lme(pixel ~ day + I(day^2), data = Pixel,
random = list(Dog = ~ day, Side = ~ 1))
summary(fm1)
plot(augPred(fm1))
I hope this helps.
Regards
Wayne
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] Behalf Of John McHenry
Sent: 13 September 2007 03:09
To: r-help@r-project.org
Subject: [R] Multivariate, multilevel regression?
Dear WizaRds,
This is mostly a statistics question, but I'm figuring that R is the right
solution (even before I start!)
I have some bio data of heart rate versus time (rats taken from resting to
maximal heart rate). I want to regress heart rate on time. The data have been
normalized such that resting heart rate is zero at time=0, so that all curves
intersect at the origin (and at the origin only). The regression function is
monotonically increasing, but little is known about it.
The data are also strictly ordinal: I have a factor such that lower-order data
*must* be on a curve that is offset beneath a higher-order curve. Don't ask
where the factors come from...but *given* these factors the assumption, or
rather the *constraint*, is that lower orders are better (lower-order rats are
fitter rats with better cardiovascular response). Most important is that these
curves do not intersect because of these factors (a fitter rat can not have a
worse heart-rate
response than a less fit rat!).
Here's a schematic to show you what I mean. It's very rough, but it gets the
point across:
seq(0, 1, len=100)
f<- 1/seq(.3,1,len=6)
windows();
plot(t, sqrt(3*t), type='n', xlab='time', ylab='heart rate');
grid()
for (i in seq(along=f))
lines(t, sqrt(f[i]*t), col=ceiling(2*i))
Now I'm wondering where I should start and I'm think that this is really not
much different from having a y_i ~ x_i | factor_i ... the i^th response curve
just like a dummy variable male/female linear regression. But in some way the
factors are related (there's a dose behind it, if ya see what I mean), so they
are not independent...they're all part of a "system" (some studies have more
"juice" overall, so the whole "system" is varying from one group of rats to the
next).
This means that in some "systems" the curves will be closer or some curves
within the system will be closer together.
Here is my question for you guys: what is the best way to model this kind of
problem, especially given that each "i^th curve" could have as few as 3 data
points? If I can only vaguely assume the type of the regression function
(monotonic, rapidly increasing from the origin kind of like sqrt(t), as above)
should it be a parametric or nonparametric regression? What about those errors?
Gee! It's hard to assume anything (there're so few of 'em and they're probably
heteroskedastic!
Where do I begin? I'd gladly accept all references, pointers, and such like.
Best, and thanks for R and R-Help, all you Core guys!!!
Jack.
-
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
