[R] Repost: Estimation when interaction is present: How do I get get the parameters from nlme?
John S. Walker
jsw9c at uic.edu
Thu Jun 15 17:36:39 CEST 2006
Gday,
This is a repost since I only had one direct reply and I remain
mystified- This
may be stupidity on my part but it may not be so simple.
In brief, my problem is I'm not sure how to extract parameter
values/effect sizes from a nonlinear
regression model with a significant interaction term.
My data sets are dose response curves (force and dose) for muscle that
also have two treatments applied
Treatment A (A- or A+) and Treatment B (B-/B+). A single muscle was
used for each experiment - a full dose response curve and one treatment
from the matrix A*B (A-/B-, A+/B-, A-/B+ and A+,B+). There are 8
replicates for each combination of treatments
We fit a dose response curve to each experiment with parameters upper,
ed50 and slope; we expect treatment A to change upper and ed50. We want
to know if treatment B blocks the effect of treatment A and if so to
what degree.
This is similar to the Ludbrook example in Venables and Ripley, however
they only had one treatment and I have two.
my approach
The dataframe is structured like this:
expt treatA treatB dose force.
1 - - 0.1 20
1 - - 0.2 40
...
4 + + 0.1 20
4 +
I used a groupedData object: mydata=groupedData(force ~ dose | expt)
I used an nlme obect to model the data as follows (pseudocode):
myfit.nlme <- nlme(force ~ ss_tpl(dose, upper, ed50,slope),
fixed=list(ed50~factor(treatA)*factor(treatB)))
The function ss_tpl is a properly debugged and fully functional
selfstarting three parameter logistic function that I wrote- no problem
here. In my analysis
I also included fixed terms for the other fit parameters; upper and
slope, but my main problem is with the
ed50 so that's all I've included here.
Running an anova on the resulting object (anova(myfit.nlme) I found the
A -/B- (control) to
be significantly different from zero, treatment A was significantly
different, treatment B had no significant
effect and there was a significant interaction between treatment A and
treatment B.
The interaction term is likely to be real. The treatments are on
sequential steps in a pathway and treatment B may be blocking the
effect of treatment A, i.e. treatment B alone has no effect because it
blocks a pathway that is not active, treatment A reduces force via this
pathway and treament B therefore blocks the effect of treatment A when
used together.
From what I understand, please correct me if I'm wrong, the parameter
estimates from summary(model.nlme) are not correct for main effects if
a significant interaction is present. For example in my data treatment
B alone has no signifcant effect in the anova but the interaction term
A:B is significant. I believe The summary estimate for B is the
estimate across all levels of A. What I want to do is pull out the
estimate for B when A is not present. I suppose I can do it manually
from the list of coefficients from nls or fit a oneway model with
treatment levels A, B, AB. But I was kind of hoping there was some
extractor function.
The reason I need this is that the co-authors want to include a table
of parameter values with std errs or confidence intervals ala:
Treat upper ed50 slope
A-/B- x x x <- shows value for comparison to control studies
A+/B- x x x <-Shows A is working0
A-/B+ x x x <- Shows B has no effect alone
A+/B + x x x <-shows B blocks A (not necessarily total)
So back to my question,How do I extract estimates of the parameters
from my model object for a
specific combination of factors including the interaction term.
i.e. what is the ed50 (and std err) for A-/B-, A+/B-, A-/B+, A+/B+ ?
I think this is a fair question and one that many biomedical scientists
would need.
More information about the R-help
mailing list