A more compact way to code factors in nls is to use the syntax factor[].
Here's an example using a simplified version of Ravi's example:
n <- 200
set.seed(123)
x <- runif(n)
a <- gl(n=2, k=n/2) # a two-level factor
eps <- rnorm(n, sd=0.5)
y <- as.numeric(a) * x^.5 + eps
nls(y ~ a[]*x^b, start=list(a=c(1,1), b=0.25))
Wed, 2009-04-15 at 10:22 -0400, Ravi Varadhan wrote:
> Hi,
>
> Here is one-way to do it (the following code shows a simulation example):
>
> n <- 200
>
> set.seed(123)
>
> x <- runif(n)
>
> f <- gl(n=2, k=n/2) # a two-level factor
>
> x1 <- x * (f == 1)
>
> x2 <- x * (f == 2)
>
> a <- c(rep(2, n/2), rep(5, n/2))
>
> b <- 0.5
>
> nsim <- 100
>
> nls.coef <- matrix(NA, nsim, 3)
>
> for (i in 1:nsim) {
>
> eps <- rnorm(n, sd=0.5)
>
> y <- a * x^b + eps
>
> ans.nls <- try(nls(y ~ a1*x1^b + a2*x2^b, start=list(a1=1, a2=1, b=0.25)),
> silent=TRUE)
>
> if (class(ans.nls) != "try-error") nls.coef[i, ] <- coef(ans.nls)
>
> }
>
> apply(nls.coef, 2, summary)
>
> Hope this helps,
> Ravi.
>
> ----------------------------------------------------------------------------
> -------
>
> Ravi Varadhan, Ph.D.
>
> Assistant Professor, The Center on Aging and Health
>
> Division of Geriatric Medicine and Gerontology
>
> Johns Hopkins University
>
> Ph: (410) 502-2619
>
> Fax: (410) 614-9625
>
> Email: rvarad...@jhmi.edu
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
> ----------------------------------------------------------------------------
> --------
>
>
> -----Original Message-----
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
> Behalf Of Manuel Gutierrez
> Sent: Wednesday, April 15, 2009 8:27 AM
> To: r-help@r-project.org
> Subject: [R] nls factor
>
> I want to fit the model y=a*x^b using nls; where "a" should be different for
> each level of a factor.
> What is the easiest way to fit it? Can i do it with nls?
> I've looked the help pages and the MASS example in page 249 but the formula
> is different and I don't know how to specify it for my model.
> Thanks,
> Manuel
>
> [[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.
--
http://mutualism.williams.edu
______________________________________________
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.
