R help - Jun 2006 - Repost: Estimation when interaction is present: How do I get get the parameters from nlme?

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.

Hi John,
I think a solution is to
1. recode A and B as a single factor, AB, with four levels,
2. define each fixed effect as a function of  AB minus the intercept  
(e.g.  ed50 ~ as.factor(AB)-1).
3. extract the tTable as a data.frame with  summary(model)$tTable.

I will be interested to see what other folks suggest.
Cheers,
Hank

  and then run the Probably not the best, nor the worst solution,  
might be to recode A and B
On Jun 15, 2006, at 11:36 AM, John S. Walker wrote:
> 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.
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting- 
> guide.html
Dr. M. Hank H. Stevens, Assistant Professor
338 Pearson Hall
Botany Department
Miami University
Oxford, OH 45056

Office: (513) 529-4206
Lab: (513) 529-4262
FAX: (513) 529-4243
http://www.cas.muohio.edu/~stevenmh/
http://www.muohio.edu/ecology/
http://www.muohio.edu/botany/
"E Pluribus Unum"