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 <-
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.466288e-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, here'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
______________________________________________
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.
