Hi Rick
Thanks to Dieter Menne I did manage to solve the problem of imposing bounds on
the parameter space duirng an nlsList fit. He suggested using optim to optimize
the parameters prior to each fit. This worked well for me as I had a customized
selfStart function that then optimized the parameters for each individual
separately.
if you rewrote your selfStart function as:
powrDpltInit <-
function(mCall, LHS, data) {
xy <- sortedXyData(mCall[["x"]],LHS,data)
min.s <- min(y)
dif.s <- max(y)-min(y)
dplt.s <- 0.5
p.s <- -.20
value <- c(min.s, dplt.s, dif.s, p.s)
#function to optimize
func1 <- function(value) {
min.s <- value[1]
dif.s <- value[2]
dplt.s <- value[3]
p.s <- value[4]
y1<-rep(0,length(xy$x)) # generate vector for predicted y (y1) to evaluate
against observed y
for(cnt in 1:length(xy$x)){
y1[cnt]<- min.s + dplt.s*x[cnt] + dif.s*x[cnt]^p.s} #predicting y1 for
values of y
evl<-sum((xy$y-y1)^2) #sum of squares is function to minimize
return(evl)}
#optimizing
oppar<-optim(c(min.s , dif.s , dplt.s ,
p.s),func1,method="L-BFGS-B",
lower=cc(0.0,0.0,0.0,-10),
control=list(maxit=2000,parscale=c(1e3,1e-5,1e3,1e-1)))
#saving optimized parameters
value<-c(oppar$par[1L],oppar$par[3L],oppar$par[2L],oppar$par[4L])
names(value) <-
mCall[c("min","dplt","dif","p")]
value
}
SSpowrDplt<-selfStart(~min + dplt*x + dif*x^p,initial=powrDpltInit,
parameters=c("min","dplt","dif","p"))
this should implement the optimization - I apologize for any typos as I was
unable to check it with appropriate data.
You may want to play with parscale so that is appropriate for your expected
parameters. I also had to increase the value for tol (tolerance) in the control
list of the nlsList call to fit for all individuals in my dataset - you may
have to play with appropriate and acceptable values of tol for yours too.
Hope this works for you
Steve
Dr. Steve Oswald
Penn State Berks
>
> Hi All.
>
> I'd like to send some control arguments to the nls function when
> performing a nlsList analysis.
>
> I'm fitting a power model to some grouped data and would like to
> impose lower bounds on the estimates using the "port" algorithm.
> Obtaining the lower bound constraint works fine with a direct call to
> nls for a single level of the grouping variable. However, the bounds
> aren't imposed when performing the nlsList analysis across the levels
> of the grouping variable. Any idea why this isn't working?
>
> # Generate example data ##########
>
> trial <- 1:100
> result <- list()
>
> for (i in 1:3) {
> min <- rnorm(max(trial),250,5)
> dif <- rnorm(max(trial),500,5)
> p <- rnorm(max(trial),-.5,.1)
> e <- rnorm(max(trial),mean=0,sd=5)
> y <- min + dif*(trial)^p + e
> result[[i]] <- data.frame(y,trial,id=i)
> }
> newdf <-do.call('rbind',result)
> df.gr <- groupedData( y ~ trial | id, data=newdf)
>
>
> ### Single unit analysis
> ........................................................................
> ### The boundary condition on the dplt parameter is enforced! ..........
>
> df.one <- subset(df.gr,id==1)
>
nls(y~SSpowrDplt(trial,min,dplt,dif,p),data=df.one,algorithm="port",lower=c(0.0,0.0,0.0,-10)>
> ...... example output.......
>Nonlinear regression model
> model: y ~ SSpowrDplt(trial, min, dplt, dif, p)
> data: df.one
> min dplt dif p
>247.052 0.000 491.965 -0.462
> residual sum-of-squares: 234322
> Algorithm "port", convergence message: relative convergence (4)
> ##################################
>
> ### nlsList analysis
>
........................................................................................
> ### Boundary condition on dplt is not enforced
> .............................................
>
> Lfit.nls &n! bsp; <-
>
nlsList(y~SSpowrDplt(trial,min,dplt,dif,p),data=df.gr,control=list(algorithm='port',lower=c(0.0,0.0,0.0,-10),maxiter=100)
>
> ...... example output.......
>Call:
> Model: y ~ SSpowrDplt(trial, min, dplt, dif, p) | id
> Data: df.gr
>
>Coefficients:
> min
> Estimate Std. Error t value Pr(>|t|)
>1 276.2354 16.16609 17.087337 1.086442e-44
>2 257.0127 20.43564 12.576694 3.390686e-30
>3 206.4017 29.01315 7.114075 7.354863e-12
> dplt
> Estimate Std. Error t value Pr(>|t|)
>1 -0.06579982 0.03848086 -1.7099365 0.0951222
>2 -0.01694362 0.04161933 -0.4071093 0.6786473
>3 0.08981518 0.04636532 1.9371199 0.0528957
> dif
> Estimate Std. Error t value Pr(>|t|)
>1 477.5049 21.89002 21.81382 ! 6.679439e-62
>2 488.7171 22.11908 22.09482 4.4662! 88e-66>3 552.7105 25.04206
> 22.07129 9.215344e-65
> p
> Estimate Std. Error t value Pr(>|t|)
>1 -0.5455936 0.06262040 -8.712713 7.615265e-16
>2 -0.4839114 0.06074282 -7.966560 1.307734e-14
>3 -0.4059903 0.05455864 -7.441355 9.297527e-13
>
>Residual standard error: 27.43384 on 888 degrees of freedom
> #####################################################
>
> If you look at the structure of Lfit.nls, it looks like the control
> arguments are passed correctly.
> str(Lfit.nls)
>
> List of 3
> $ 1:List of 6
> ..$ control :List of 7
> .. ..$ maxiter : num 100
> .. ..$ tol : num 1e-05
> .. ..$ minFactor: num 0.000977
> .. ..$ printEval: logi FALSE
> .. ..$ warnOnly : logi FALSE
> .. ..$ algorithm: chr "port"
> .. ..$ lower : num [1:4] 0 0 0 -10
>
>
>
> If it helps, he! re's the selfStart function that I'm using....
> powrDpltInit <-
> function(mCall, LHS, data) {
> xy <- sortedXyData(mCall[["x"]],LHS,data)
> min.s <- min(y)
> dif.s <- max(y)-min(y)
> dplt.s <- 0.5
> p.s <- -.20
> value <- c(min.s, dplt.s, dif.s, p.s)
> names(value) <-
mCall[c("min","dplt","dif","p")]
> value
> }
>
> SSpowrDplt<-selfStart(~min + dplt*x + dif*x^p,initial=powrDpltInit,
>
parameters=c("min","dplt","dif","p")
>
>
>
> Thanks for your help!
>
> Rick DeShon
>
> ______________________________________________
>
[[alternative HTML version deleted]]
