I am interested in testing two similar nls models to determine if the
lines are statistically different when fitted with two different data
sets; one corn, another soybean. I know I can do this in linear models
by testing for interactions. See Introductory Statistics with R by
Dallgaard p212-218 for an example. I have two different data sets I am
comparing to lai. ci.re should have very little difference between corn
and soybean, while ci.gr, there should be a difference. If I use the
simplistic form described in Dallgaard (see example in code below) it
does not work correctly. What do I need to do to test for this
interaction? Thank you!
My data is located here: ftp://snrftp.unl.edu/Outgoing/example/
Here is my example code:
###load data###
data <- read.csv("example.csv", header=TRUE)
eq <- function(x,a,b,c) {a+b*exp(x^c)}
##create non-linear models##
nls.ci.re <- nls(lai~eq(ci.re,a,b,c),data=data, start=c(a=1,b=1,c=1))
nls.ci.gr <- nls(lai~eq(ci.gr,a,b,c),data=data, start=c(a=1,b=1,c=1))
#create non-linear models for corn#
nls.ci.re.corn <-
nls(lai~eq(ci.re,a,b,c),data=subset(data,crop=="corn"),
start=c(a=1,b=1,c=1))
nls.ci.gr.corn <-
nls(lai~eq(ci.gr,a,b,c),data=subset(data,crop=="corn"),
start=c(a=1,b=1,c=1))
#create non-linear models for soybean#
nls.ci.re.soybean <-
nls(lai~eq(ci.re,a,b,c),data=subset(data,crop=="soybean"),
start=c(a=1,b=1,c=1))
nls.ci.gr.soybean <-
nls(lai~eq(ci.gr,a,b,c),data=subset(data,crop=="soybean"),
start=c(a=1,b=1,c=1))
#test interactions according to Introductory Statistics with R by
Dalgaard p. 213#
corn.soybean.interactions.ci.re.a <-
abs((summary(nls.ci.re.corn)$coeff[1,1]-
summary(nls.ci.re.soybean)$coeff[1,1])/
sqrt(summary(nls.ci.re.corn)$coeff[1,2]^2+
summary(nls.ci.re.soybean)$coeff[1,2]^2))
corn.soybean.interactions.ci.gr.a <-
abs((summary(nls.ci.gr.corn)$coeff[1,1]-
summary(nls.ci.gr.soybean)$coeff[1,1])/
sqrt(summary(nls.ci.gr.corn)$coeff[1,2]^2+
summary(nls.ci.gr.soybean)$coeff[1,2]^2))
corn.soybean.interactions.ci.re.a.p.value <-
pt(corn.soybean.interactions.ci.re.a,df=summary(nls.ci.re.corn)$df[2],lower.tail=FALSE)
corn.soybean.interactions.ci.gr.a.p.value <-
pt(corn.soybean.interactions.ci.gr.a,df=summary(nls.ci.gr.corn)$df[2],lower.tail=FALSE)
Nonlinear models are an entirely different kettle of fish then linear models. You need to specify exactly how the different crops affect the parameters in your growth model. I suggest you consult a local statistician for help (this is not an R question). Incidentally, depending on the nature of your data -- does it consist of repeated measures of individual plants, for example? -- you may need nlme or nlmer (or other special packages) to fit the data and do inference -- nls makes independence assumptions that aren't appropriate for such data. Your local statistician can help here, too, if needed. But be warned, inference is inherently difficult in such situations. -- Bert On Tue, Feb 8, 2011 at 1:12 PM, Anthony Lawrence Nguy-Robertson <anthony.robertson at huskers.unl.edu> wrote:> I am interested in testing two similar nls models to determine if the lines > are statistically different when fitted with two different data sets; one > corn, another soybean. I know I can do this in linear models by testing for > interactions. See Introductory Statistics with R by Dallgaard p212-218 for > an example. I have two different data sets I am comparing to lai. ci.re > should have very little difference between corn and soybean, while ci.gr, > there should be a difference. If I use the simplistic form described in > Dallgaard (see example in code below) it does not work correctly. What do I > need to do to test for this interaction? Thank you! > > My data is located here: ftp://snrftp.unl.edu/Outgoing/example/ > > Here is my example code: > > ###load data### > data <- read.csv("example.csv", header=TRUE) > eq <- function(x,a,b,c) {a+b*exp(x^c)} > > ##create non-linear models## > nls.ci.re <- nls(lai~eq(ci.re,a,b,c),data=data, start=c(a=1,b=1,c=1)) > nls.ci.gr <- nls(lai~eq(ci.gr,a,b,c),data=data, start=c(a=1,b=1,c=1)) > > #create non-linear models for corn# > nls.ci.re.corn <- nls(lai~eq(ci.re,a,b,c),data=subset(data,crop=="corn"), > start=c(a=1,b=1,c=1)) > nls.ci.gr.corn <- nls(lai~eq(ci.gr,a,b,c),data=subset(data,crop=="corn"), > start=c(a=1,b=1,c=1)) > > #create non-linear models for soybean# > nls.ci.re.soybean <- > nls(lai~eq(ci.re,a,b,c),data=subset(data,crop=="soybean"), > start=c(a=1,b=1,c=1)) > nls.ci.gr.soybean <- > nls(lai~eq(ci.gr,a,b,c),data=subset(data,crop=="soybean"), > start=c(a=1,b=1,c=1)) > > #test interactions according to Introductory Statistics with R by Dalgaard > p. 213# > corn.soybean.interactions.ci.re.a <- > ? ?abs((summary(nls.ci.re.corn)$coeff[1,1]- > ? ?summary(nls.ci.re.soybean)$coeff[1,1])/ > ? ?sqrt(summary(nls.ci.re.corn)$coeff[1,2]^2+ > ? ?summary(nls.ci.re.soybean)$coeff[1,2]^2)) > > corn.soybean.interactions.ci.gr.a <- > ? ?abs((summary(nls.ci.gr.corn)$coeff[1,1]- > ? ?summary(nls.ci.gr.soybean)$coeff[1,1])/ > ? ?sqrt(summary(nls.ci.gr.corn)$coeff[1,2]^2+ > ? ?summary(nls.ci.gr.soybean)$coeff[1,2]^2)) > > corn.soybean.interactions.ci.re.a.p.value <- > pt(corn.soybean.interactions.ci.re.a,df=summary(nls.ci.re.corn)$df[2],lower.tail=FALSE) > corn.soybean.interactions.ci.gr.a.p.value <- > pt(corn.soybean.interactions.ci.gr.a,df=summary(nls.ci.gr.corn)$df[2],lower.tail=FALSE) > > ______________________________________________ > R-help at 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. >-- Bert Gunter Genentech Nonclinical Biostatistics
On Tue, Feb 8, 2011 at 4:12 PM, Anthony Lawrence Nguy-Robertson <anthony.robertson at huskers.unl.edu> wrote:> I am interested in testing two similar nls models to determine if the lines > are statistically different when fitted with two different data sets; one > corn, another soybean. I know I can do this in linear models by testing for > interactions. See Introductory Statistics with R by Dallgaard p212-218 for > an example. I have two different data sets I am comparing to lai. ci.re > should have very little difference between corn and soybean, while ci.gr, > there should be a difference. If I use the simplistic form described in > Dallgaard (see example in code below) it does not work correctly. What do I > need to do to test for this interaction? Thank you! > > My data is located here: ftp://snrftp.unl.edu/Outgoing/example/ > > Here is my example code: > > ###load data### > data <- read.csv("example.csv", header=TRUE) > eq <- function(x,a,b,c) {a+b*exp(x^c)} > > ##create non-linear models## > nls.ci.re <- nls(lai~eq(ci.re,a,b,c),data=data, start=c(a=1,b=1,c=1)) > nls.ci.gr <- nls(lai~eq(ci.gr,a,b,c),data=data, start=c(a=1,b=1,c=1)) > > #create non-linear models for corn# > nls.ci.re.corn <- nls(lai~eq(ci.re,a,b,c),data=subset(data,crop=="corn"), > start=c(a=1,b=1,c=1)) > nls.ci.gr.corn <- nls(lai~eq(ci.gr,a,b,c),data=subset(data,crop=="corn"), > start=c(a=1,b=1,c=1)) > > #create non-linear models for soybean# > nls.ci.re.soybean <- > nls(lai~eq(ci.re,a,b,c),data=subset(data,crop=="soybean"), > start=c(a=1,b=1,c=1)) > nls.ci.gr.soybean <- > nls(lai~eq(ci.gr,a,b,c),data=subset(data,crop=="soybean"), > start=c(a=1,b=1,c=1)) > > #test interactions according to Introductory Statistics with R by Dalgaard > p. 213# > corn.soybean.interactions.ci.re.a <- > ? ?abs((summary(nls.ci.re.corn)$coeff[1,1]- > ? ?summary(nls.ci.re.soybean)$coeff[1,1])/ > ? ?sqrt(summary(nls.ci.re.corn)$coeff[1,2]^2+ > ? ?summary(nls.ci.re.soybean)$coeff[1,2]^2)) > > corn.soybean.interactions.ci.gr.a <- > ? ?abs((summary(nls.ci.gr.corn)$coeff[1,1]- > ? ?summary(nls.ci.gr.soybean)$coeff[1,1])/ > ? ?sqrt(summary(nls.ci.gr.corn)$coeff[1,2]^2+ > ? ?summary(nls.ci.gr.soybean)$coeff[1,2]^2)) > > corn.soybean.interactions.ci.re.a.p.value <- > pt(corn.soybean.interactions.ci.re.a,df=summary(nls.ci.re.corn)$df[2],lower.tail=FALSE) > corn.soybean.interactions.ci.gr.a.p.value <- > pt(corn.soybean.interactions.ci.gr.a,df=summary(nls.ci.gr.corn)$df[2],lower.tail=FALSE)There is an anova.nls method: anova(model1, model2) -- Statistics & Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at gmail.com
On Tue, Feb 8, 2011 at 11:43 PM, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:> On Tue, Feb 8, 2011 at 4:12 PM, Anthony Lawrence Nguy-Robertson > <anthony.robertson at huskers.unl.edu> wrote: >> I am interested in testing two similar nls models to determine if the lines >> are statistically different when fitted with two different data sets; one >> corn, another soybean. I know I can do this in linear models by testing for >> interactions. See Introductory Statistics with R by Dallgaard p212-218 for >> an example. I have two different data sets I am comparing to lai. ci.re >> should have very little difference between corn and soybean, while ci.gr, >> there should be a difference. If I use the simplistic form described in >> Dallgaard (see example in code below) it does not work correctly. What do I >> need to do to test for this interaction? Thank you! >> >> My data is located here: ftp://snrftp.unl.edu/Outgoing/example/ >> >> Here is my example code: >> >> ###load data### >> data <- read.csv("example.csv", header=TRUE) >> eq <- function(x,a,b,c) {a+b*exp(x^c)} >> >> ##create non-linear models## >> nls.ci.re <- nls(lai~eq(ci.re,a,b,c),data=data, start=c(a=1,b=1,c=1)) >> nls.ci.gr <- nls(lai~eq(ci.gr,a,b,c),data=data, start=c(a=1,b=1,c=1)) >> >> #create non-linear models for corn# >> nls.ci.re.corn <- nls(lai~eq(ci.re,a,b,c),data=subset(data,crop=="corn"), >> start=c(a=1,b=1,c=1)) >> nls.ci.gr.corn <- nls(lai~eq(ci.gr,a,b,c),data=subset(data,crop=="corn"), >> start=c(a=1,b=1,c=1)) >> >> #create non-linear models for soybean# >> nls.ci.re.soybean <- >> nls(lai~eq(ci.re,a,b,c),data=subset(data,crop=="soybean"), >> start=c(a=1,b=1,c=1)) >> nls.ci.gr.soybean <- >> nls(lai~eq(ci.gr,a,b,c),data=subset(data,crop=="soybean"), >> start=c(a=1,b=1,c=1)) >> >> #test interactions according to Introductory Statistics with R by Dalgaard >> p. 213# >> corn.soybean.interactions.ci.re.a <- >> ? ?abs((summary(nls.ci.re.corn)$coeff[1,1]- >> ? ?summary(nls.ci.re.soybean)$coeff[1,1])/ >> ? ?sqrt(summary(nls.ci.re.corn)$coeff[1,2]^2+ >> ? ?summary(nls.ci.re.soybean)$coeff[1,2]^2)) >> >> corn.soybean.interactions.ci.gr.a <- >> ? ?abs((summary(nls.ci.gr.corn)$coeff[1,1]- >> ? ?summary(nls.ci.gr.soybean)$coeff[1,1])/ >> ? ?sqrt(summary(nls.ci.gr.corn)$coeff[1,2]^2+ >> ? ?summary(nls.ci.gr.soybean)$coeff[1,2]^2)) >> >> corn.soybean.interactions.ci.re.a.p.value <- >> pt(corn.soybean.interactions.ci.re.a,df=summary(nls.ci.re.corn)$df[2],lower.tail=FALSE) >> corn.soybean.interactions.ci.gr.a.p.value <- >> pt(corn.soybean.interactions.ci.gr.a,df=summary(nls.ci.gr.corn)$df[2],lower.tail=FALSE) > > There is an anova.nls method: anova(model1, model2)Just to be clear the two models would each include both groups -- one model would assume the parameters are the same for the two groups so it would have 3 parameters and the other model would allow them to be different (up to 6 parameters depending on how many parameters you wish to allow to be different between the two groups) -- these are not the models shown above. -- Statistics & Software Consulting GKX Group, GKX Associates Inc. tel: 1-877-GKX-GROUP email: ggrothendieck at gmail.com
Thank you R-forum for you generous help. Gabor Grothendieck, I am not sure if anova in the form that you suggested is the most appropriate (This is probably more of a statistics related, rather than R related at this point). The way I understand anova is that you are testing the variance between the models. I know that the variance in the 'corn' models is greater than 'soybean' due to a biological reason. I don't think this approach is correct in my case since the anova analysis is comparing variance between different data sets. However, using the approach suggested by Derek Ogle seems to produce results that make more sense since the anova analysis is between the entire data set and the influence of crop type on that model. For the couple of models I tested thus far, I think this is the approach I need to take. I am not a statitician, so I hope my descriptions make sense. Again I appreciate everyone's help in this matter. Thanks, Tony On 2/9/2011 7:53 AM, Derek Ogle wrote:> Anthony, > > The link below shows how I did a similar analysis but for a different nls model. Perhaps it will be useful to you ... > > http://www.ncfaculty.net/dogle/fishR/gnrlex/VonBertalanffyCompare/VonBertalanffyCompare.pdf > > Good luck. > > >> -----Original Message----- >> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r- >> project.org] On Behalf Of Gabor Grothendieck >> Sent: Tuesday, February 08, 2011 10:43 PM >> To: Anthony Lawrence Nguy-Robertson >> Cc: r-help at r-project.org >> Subject: Re: [R] Interactions in a nls model >> >> On Tue, Feb 8, 2011 at 4:12 PM, Anthony Lawrence Nguy-Robertson >> <anthony.robertson at huskers.unl.edu> wrote: >> >>> I am interested in testing two similar nls models to determine if the >>> >> lines >> >>> are statistically different when fitted with two different data sets; >>> >> one >> >>> corn, another soybean. I know I can do this in linear models by >>> >> testing for >> >>> interactions. See Introductory Statistics with R by Dallgaard p212- >>> >> 218 for >> >>> an example. I have two different data sets I am comparing to lai. >>> >> ci.re >> >>> should have very little difference between corn and soybean, while >>> >> ci.gr, >> >>> there should be a difference. If I use the simplistic form described >>> >> in >> >>> Dallgaard (see example in code below) it does not work correctly. >>> >> What do I >> >>> need to do to test for this interaction? Thank you! >>> >>> My data is located here: ftp://snrftp.unl.edu/Outgoing/example/ >>> >>> Here is my example code: >>> >>> ###load data### >>> data<- read.csv("example.csv", header=TRUE) >>> eq<- function(x,a,b,c) {a+b*exp(x^c)} >>> >>> ##create non-linear models## >>> nls.ci.re<- nls(lai~eq(ci.re,a,b,c),data=data, start=c(a=1,b=1,c=1)) >>> nls.ci.gr<- nls(lai~eq(ci.gr,a,b,c),data=data, start=c(a=1,b=1,c=1)) >>> >>> #create non-linear models for corn# >>> nls.ci.re.corn<- >>> >> nls(lai~eq(ci.re,a,b,c),data=subset(data,crop=="corn"), >> >>> start=c(a=1,b=1,c=1)) >>> nls.ci.gr.corn<- >>> >> nls(lai~eq(ci.gr,a,b,c),data=subset(data,crop=="corn"), >> >>> start=c(a=1,b=1,c=1)) >>> >>> #create non-linear models for soybean# >>> nls.ci.re.soybean<- >>> nls(lai~eq(ci.re,a,b,c),data=subset(data,crop=="soybean"), >>> start=c(a=1,b=1,c=1)) >>> nls.ci.gr.soybean<- >>> nls(lai~eq(ci.gr,a,b,c),data=subset(data,crop=="soybean"), >>> start=c(a=1,b=1,c=1)) >>> >>> #test interactions according to Introductory Statistics with R by >>> >> Dalgaard >> >>> p. 213# >>> corn.soybean.interactions.ci.re.a<- >>> abs((summary(nls.ci.re.corn)$coeff[1,1]- >>> summary(nls.ci.re.soybean)$coeff[1,1])/ >>> sqrt(summary(nls.ci.re.corn)$coeff[1,2]^2+ >>> summary(nls.ci.re.soybean)$coeff[1,2]^2)) >>> >>> corn.soybean.interactions.ci.gr.a<- >>> abs((summary(nls.ci.gr.corn)$coeff[1,1]- >>> summary(nls.ci.gr.soybean)$coeff[1,1])/ >>> sqrt(summary(nls.ci.gr.corn)$coeff[1,2]^2+ >>> summary(nls.ci.gr.soybean)$coeff[1,2]^2)) >>> >>> corn.soybean.interactions.ci.re.a.p.value<- >>> >>> >> pt(corn.soybean.interactions.ci.re.a,df=summary(nls.ci.re.corn)$df[2],l >> ower.tail=FALSE) >> >>> corn.soybean.interactions.ci.gr.a.p.value<- >>> >>> >> pt(corn.soybean.interactions.ci.gr.a,df=summary(nls.ci.gr.corn)$df[2],l >> ower.tail=FALSE) >> >> There is an anova.nls method: anova(model1, model2) >> >> >> -- >> Statistics& Software Consulting >> GKX Group, GKX Associates Inc. >> tel: 1-877-GKX-GROUP >> email: ggrothendieck at gmail.com >> >> ______________________________________________ >> R-help at 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. >> > . > >-- Anthony L. Nguy-Robertson Doctoral Student University of Nebraska-Lincoln School of Natural Resources Center for Advanced Land Management Information Technologies (CALMIT) 223 Hardin Hall 3310 Holdrege Street Lincoln, NE 68583 Office: 402-472-2565 Fax: 402-472-2946 Email: anguy-robertson at unlnotes.unl.edu