[R] nls factor

[Manuel Morales]

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: rvaradhan at jhmi.edu
> 
> Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
> 
>  
> 
> ----------------------------------------------------------------------------
> --------
> 
> 
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
> Behalf Of Manuel Gutierrez
> Sent: Wednesday, April 15, 2009 8:27 AM
> To: r-help at 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 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.
> 
> ______________________________________________
> 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.
-- 
http://mutualism.williams.edu