[R] get optim results into a model object
Ravi Varadhan
RVaradhan at jhmi.edu
Wed Apr 8 00:17:49 CEST 2009
Gene,
I agree that by using "nls" you are able to also enforce box-constraints on
the betas. There is one more issue. My trick doesn't really force the sum
to be 1. It is very close, but not exactly so. You can use "weights" to
achieve this. Here you go:
# analysis without "weights"
ans2 <- nls(aformula,data=df,,algorithm='port',
start=list(b0=0,b1=1/3,b2=1/3,b3=1/3),
lower=c(-Inf,rep(0,3)),upper=c(Inf,rep(1,3)), trace=T)
> sum(summary(ans2)$coef[2:4,1])
[1] 0.9999997
>
# Now, let us use a large weight for the first "artificial" observation:
wts <- c(100, rep(1, length(y)))
ans2 <- nls(aformula,data=df,,algorithm='port',
start=list(b0=0,b1=1/3,b2=1/3,b3=1/3),
lower=c(-Inf,rep(0,3)),upper=c(Inf,rep(1,3)), weights=wts, trace=T)
> sum(summary(ans2)$coef[2:4,1])
[1] 1
>
I know the difference is insignificant, in this example, but in general,
using a large weight for the first equation is a good idea.
Best,
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 Gene Leynes
Sent: Tuesday, April 07, 2009 5:46 PM
To: Ravi Varadhan
Cc: r-help at r-project.org
Subject: Re: [R] get optim results into a model object
Very nice trick!!! Thank you!
When you combine your trick with the "nls" function, you get the whole
thing:
#(I have also simplified the model since the first post)
###################################################################
#Original Method: uses optim, but doesn't create "model":
###################################################################
x=matrix(c(0.04, 0.08, 0.01, 0.0398, 0.081, 0.00915, 0.04057, 0.07944,
0.00994, 0.04137, 0.07949, 0.01132, 0.0421, 0.08273, 0.00947, 0.04237,
0.08058, 0.00969, 0.0425, 0.07952, 0.00919, 0.04305, 0.07717, 0.00908,
0.04319, 0.07576, 0.0061, 0.04298, 0.07557, 0.00539, 0.04287, 0.07244,
0.00542, 0.04318, 0.071, 0.00388, 0.04348, 0.07053, 0.00375, 0.04335,
0.07075, 0.00364, 0.04386, 0.07481, 0.00296),ncol=3,byrow=T)
x=cbind(1,x)
y=matrix(c(0.02847,0.02831,0.03255,0.02736,0.03289,0.02774,0.03192,0.03141,
0.02927,0.02648,0.02807,0.03047,0.03046,0.02541,0.02729))
plot(x=1:15,ylim=c(max(x[,-1]),min(x[,-1])))
matlines(cbind(x[,-1],y))
# "normalization" fun to make the rest of the x's add up to 1
normalize=function(x)c(x[1], x[2:length(x)]/sum(x[2:length(x)]))
# objective function:
fn=function(ws){
ws=normalize(ws)
sum((x %*% ws - yhat)^2)}
llim = c(-Inf,rep(0,ncol(x)-1)) # alpha is -Inf to Inf ulim = c(
Inf,rep(1,ncol(x)-1)) # betas are 0 to 1
th=matrix(c(0,rep(1/3,3)))
o = optim(th, fn, method="L-BFGS-B",lower=llim, upper=ulim) NormPar =
normalize(o$par) cat('pars:', NormPar,'\n', 'tot:', sum(NormPar), '\n')
###################################################################
# Dr. Varadhan's Method to make betas sum to 1:
###################################################################
xnew=rbind(c(0, 1, 1, 1),x)
ynew=rbind(1,y)
ans <- lm(ynew ~ xnew - 1)
summary(ans)
sum(ans$resid^2) # compare with the objective function value from
cat('pars:', ans$coef,'\n', 'tot:', sum(ans$coef[2:4]), '\n')
###################################################################
# Dr. Varadhan's Method WITH nls Method
###################################################################
xnew=rbind(c(0, 1, 1, 1),x)
ynew=rbind(1,y)
df=data.frame(cbind(ynew,xnew))
colnames(df)=c('y','x0','x1','x2','x3')
aformula = y ~ b0*x0 + b1*x1 + b2*x2 + b3*x3
ans2 <- nls(aformula,data=df,,algorithm='port',
start=list(b0=0,b1=1/3,b2=1/3,b3=1/3),
lower=c(-Inf,rep(0,3)),upper=c(Inf,rep(1,3)),trace=T)
summary(ans2)
sum(residuals(ans2)^2) # compare with the objective function value from
cat('pars:', coef(ans2),'\n', 'tot:', sum(coef(ans2)[2:4]), '\n')
On Tue, Apr 7, 2009 at 12:37 PM, Ravi Varadhan <RVaradhan at jhmi.edu> wrote:
> Hi Gene,
>
> Try the following approach using lm():
>
> set.seed(121)
> x1=.04
> for (i in 1:14) x1[i+1]=x1[i]*(1+rnorm(1)*.008)+.00025
> x2=.08
> for (i in 1:14) x2[i+1]=x2[i]*(1+rnorm(1)*.03)-.0018
> x3=.01
> for (i in 1:14) x3[i+1]=x3[i]*(1+rnorm(1)*.15)-.0008
>
> b=matrix(c(0.6,0.0,0.4)) # why don't you assume a value for intercept?
> x=matrix(cbind(x1,x2,x3),ncol=3)
>
> y=x%*%b # the 'real' y
> yhat=y+runif(15)*.006 # the observed y
> x=cbind(1,x)
>
> # here is the simple trick: add a row to X matrix and an element to
> yhat vector to impose constraints # xadd <- c(0, 1, 1, 1) xnew <-
> rbind(xadd, x) ynew <- c(1, yhat)
>
> ans <- lm(ynew ~ xnew - 1)
> summary(ans)
> sum(ans$coef[2:4])
> sum(ans$resid^2) # compare with the objective function value
> from your approach
> > sum(ans$resid^2)
> [1] 2.677474e-05
> >
>
> > o$value
> [1] 2.718646e-05
> >
>
> Note that the sum of squared residuals from lm() is smaller than your
> value from optim(). Although, this approach works well in your
> example, it does not guarantee that the coefficients are between 0 and 1.
>
> 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 Gene Leynes
> Sent: Tuesday, April 07, 2009 12:17 PM
> To: r-help at r-project.org
> Subject: [R] get optim results into a model object
>
> Hello all, I have an optimization routine that is giving me good
> results, but the results are not in the nice "model" format like "lm".
> How can I get optim results into a model so that I can use the clever
> 'fitted', 'residuals', and 'summary' functions?
>
> Using optim is the only way that I was able to make a model that
> 1) sums the betas to 1,
> 2) constrains the betas to positive numbers less than 1
> 3) does not constrain alpha
> (The constrOptim cousin wasn't very accurate, and was very slow.)
>
> Here is an example of some code, the results of which I would like to
> get into a model with the form y ~ alpha + REALPAR * x where 'REALPAR'
> is the "normalized" output at the very end
>
> many thanks!!!!
> ++++++++++++++++Code Below++++++++++++++++++++++++++++
>
> set.seed(121)
> x1=.04
> for (i in 1:14) x1[i+1]=x1[i]*(1+rnorm(1)*.008)+.00025
> x2=.08
> for (i in 1:14) x2[i+1]=x2[i]*(1+rnorm(1)*.03)-.0018
> x3=.01
> for (i in 1:14) x3[i+1]=x3[i]*(1+rnorm(1)*.15)-.0008
>
> b=matrix(c(0.6,0.0,0.4))
> x=matrix(cbind(x1,x2,x3),ncol=3)
> y=x%*%b # the 'real' y
> yhat=y+runif(15)*.006 # the observed y
>
> plot(x=1:15,ylim=c(min(x1,x2,x3),max(x1,x2,x3)))
> matlines(cbind(x,y,yhat))
>
> # Add a constant to x (for alpha)
> x=cbind(1,x)
> # "normalization" fun to make the rest of the x's add up to 1
> normalize=function(x)c(x[1], x[2:length(x)]/sum(x[2:length(x)]))
> # objective function:
> fn=function(ws){
> ws=normalize(ws)
> sum((x %*% ws - yhat)^2)}
>
> llim = c(-Inf,rep(0,ncol(x)-1)) # alpha (col 1 of 'x') is -Inf to Inf
> ulim = c( Inf,rep(1,ncol(x)-1)) # betas (cols 2:4 of 'x') are 0 to 1
> th=matrix(c(0,rep(1/3,3)))
> o = optim(th, fn, method="L-BFGS-B",lower=llim, upper=ulim) o REALPAR
> =
> normalize(o$par) REALPAR
>
> [[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.
>
>
[[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.
More information about the R-help
mailing list