5
$\begingroup$

Problem

I would like to estimate $\lambda$ in the fitted line $x^\lambda$, where $x \in [0, 1]$.

Note that the following R code generates "concave" growth as x increases from $0$ to $1$.

Lambda = 1/2.42

x = rbeta(1e4, shape1=2,shape2=2)
y = x^Lambda + rnorm(1e4, sd=.1)

plot(x,y)

Try

My approach is to build a loss function, the L2 loss to find the optimal $\lambda$, like the following R code.

SumSq = function(lam) sum((y - x^lam )^2)
optimize(SumSq, c(0,1), tol=1e-4, maximum = F )

This code assumes I know $\lambda \in (0,1)$, and optimize gives me a quite accurate result.

But this approach cannot give any statistical asymptotic results such as CI, which is useful information if I would like to consider the uncertainty.

Is there any standard way to do it?

$\endgroup$
  • $\begingroup$ i really dont know why someone would down vote and not leave a comment. This is a reasonable question with a minimally reproducible example. $\endgroup$ – forecaster Mar 7 at 0:48
  • 2
    $\begingroup$ @forecaster I agree and found your answer to this question to be edifying, so I am similarly puzzled. $\endgroup$ – James Phillips Mar 7 at 2:52
  • 3
    $\begingroup$ The natural approach to estimations with this problem would involve calling a nonlinear least squares routine (nls in R - there are a number of somewhat more robust alternatves, but this is usually sufficient). $\endgroup$ – Glen_b -Reinstate Monica Mar 7 at 5:25
  • 3
    $\begingroup$ When the model is expressed in the form $Y= \exp(\lambda \log x) + \varepsilon,$ this question is seen to be the same as many other questions asked (and answered) on this site about nonlinear least-squares fitting: this gives you an additional set of resources to consult. $\endgroup$ – whuber Mar 7 at 16:29
  • $\begingroup$ "But this approach cannot give any statistical asymptotic results such as CI, which is useful information if I would like to consider the uncertainty." Is this the core of the question? This is not so clear from the title which seems to point to point estimates (pun intended). $\endgroup$ – Sextus Empiricus Mar 7 at 19:43
7
$\begingroup$

Firstly I would suggest you read an excellent book by Ben Bolker entitled Ecological Models and Data in R. Written for ecologists , I think this is one of the best books on practical data analysis regardless of background, I'm an engineer and I have used it a lot. I'm sure there are other mathematical statistics book, however this book is by far the most practical book that I have read. He also has a package BBMLE which you might want to check. The book unlike any other goes in to maximum likelihood estimation, profile likelihood and confidence interval estimation and all that.

Coming back to your problem, you need to write a log-likelihood of the function that you are trying to fit the data. Obviously you are assuming the error to be normally distributed as in the simulation, so the Log likelihood is:

$Log\ Likelihood (Lambda,sd) = -\frac{n}{2} log(sd) - \frac{n}{2} log(2\pi)- \frac{(y-x^{Lambda})^2}{2sd} $

You need to maximize the above equation using an optimization routine such as optim in R. Use the Hessian matrix from the optimization to assess the uncertainty of your estimated log likelihood function i.e., how steep or how flat the curvature of your function is at the optimal point. If the function is steep which implies less uncertainties at optimal point you would have a tighter confidence band on your parameter estimates, on the other hand if its flat you would have a wider confidence interval. Hessian, Fisher Information matrix would help you calculate the standard error and confidence interval.

Here is how you do it in R:

    set.seed(8345)

Lambda = 1/2.42

x = rbeta(1e4, shape1=2,shape2=2)
y = x^Lambda + rnorm(1e4,mean = 0, sd=.1)

plot(x,y)

## Write Log Likelihood function

log.lik <- function(theta,y,x){

  Lam <- theta[1]
  sigma2 <- theta[2]

  # sample size
  n <-  length(y)

  #error
  e<-y-(x^Lam)

  #log likelihood
  logl<- -.5*n*log(2*pi)-.5*n*log(sigma2)-((t(e)%*%e)/(2*sigma2))

  return(-logl) # R optim does minimize so to maximize  multiply by -1
}

## Estimate Paramters thru maximum likelihood

max.lik <- optim(c(1,1), fn=log.lik, method = "L-BFGS-B", lower = c(0.00001,0.00001), hessian = T,y=y,x=x)

# Lambda
Lam <- max.lik$par[1]
#0.4107119

#Fisher Information MAtrix
fisher_info<-solve(max.lik$hessian)
prop_sigma<-sqrt(diag(fisher_info))

## Estimate 95% Confidence Interval
upper<-max.lik$par+1.96*prop_sigma
lower<-max.lik$par-1.96*prop_sigma


interval<-data.frame(Parameter = c("Lambda","sd"),value=max.lik$par, lower=lower, upper=upper)
interval
$\endgroup$
  • 3
    $\begingroup$ This is a generalized linear model so it seems more obvious to use the glm function and avoid a general non-linear optimizer with numerical differentiation. $\endgroup$ – Benjamin Christoffersen Mar 7 at 21:10
5
$\begingroup$

Is there any standard way to do it?

If you think that following is a good approximation (the model you simulate from as an example)

$$y_i = x_i^\lambda +\epsilon_i, \qquad \epsilon_i\sim N(0,\sigma^2)$$

i.e.,

$$\log(E(y_i))=\lambda \log x_i$$

then you can use glm with family = gaussian("log") which is exactly this model

lambda <- 1/2.42

set.seed(49564503)
n <- 1e4
x <-  rbeta(n, shape1 = 2,shape2 = 2)
y <-  x^lambda + rnorm(n, sd = .1)

fit <- glm(y ~ log(x) - 1, family = gaussian("log"), start =  c(0, 1))
summary(fit)
#R 
#R Call:
#R glm(formula = y ~ log(x) - 1, family = gaussian("log"), start = 1)
#R 
#R Deviance Residuals: 
#R      Min        1Q    Median        3Q       Max  
#R -0.42452  -0.06767   0.00090   0.07020   0.34418  
#R 
#R Coefficients:
#R        Estimate Std. Error t value Pr(>|t|)    
#R log(x) 0.410666   0.001807   227.3   <2e-16 ***
#R ---
#R Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#R 
#R (Dispersion parameter for gaussian family taken to be 0.01033428)
#R 
#R     Null deviance: 1057.81  on 10000  degrees of freedom
#R Residual deviance:  103.33  on  9999  degrees of freedom
#R AIC: -17341
#R 
#R Number of Fisher Scoring iterations: 5

Notice that both the dispersion parameter and coefficient estimate match as expected. Now, you can make confidence interval with the standard error above or by using confint.

$\endgroup$
  • 1
    $\begingroup$ +1 There seem to be a - 1 missing in the model formula although it's there in the output from summary(fit). $\endgroup$ – Jarle Tufto Mar 7 at 21:51
  • $\begingroup$ Thanks. I forgot to add that. I would prefer to keep the intercept but then the model would be more general than the model the OP simulates from. $\endgroup$ – Benjamin Christoffersen Mar 7 at 22:17
  • 1
    $\begingroup$ Your analysis confuses the model $E[\log y] = \lambda \log x$ with the different model $\log(E[y]) = \lambda \log x.$ When the variance of $\epsilon$ is much smaller than $x^{2\lambda}$ it's probably an OK approximation, but otherwise it's important to pay attention to the distinction--and to offer a different solution. $\endgroup$ – whuber Mar 8 at 0:28
  • 1
    $\begingroup$ @whuber doesn't the log link function express the following? $$\log(E[y]) = X\beta$$ Then using $X = x^\prime = \log(x)$ you get: $$\log(E[y]) = \beta x^\prime \quad \rightarrow \quad E[y] = e^{ \beta x^\prime} = e^{ \beta \log(x)} = x^\beta $$ $\endgroup$ – Sextus Empiricus Mar 8 at 0:41
  • 2
    $\begingroup$ @whuber, I do not see why $E [\log (y)] $ matters. The OP is using the model $$y|x \sim N (x^\lambda, \sigma^2) $$ where $\mu = x^\lambda = \exp ( \lambda \log (x))$. So $$\log (E [y|x]) = \log (\mu) = \lambda \log (x) $$ $\endgroup$ – Sextus Empiricus Mar 8 at 7:23
1
$\begingroup$

If you know how y = f(x) you can try to linearize f(x) and then use a regression. In this case, take logs of both sides such that:

log(y) = lambda log(x).

Lambda = 1/2.42
x = rbeta(1e4, shape1=2,shape2=2)
y = x^Lambda + rnorm(1e4, sd=.1)

plot(x,y)
reg <- lm(log(y)~ log(x))
summary(reg)
Lambda

You'll note that that your estimate of lambda has confidence intervals and contains the true lambda.

Note: This example has low/well behaved error such that taking the log of the error doesnt change the error that much, log(N(0,small sigma)) ~ N(0,sigma). If the error is sufficiently large enough/wonky enough, you can try to use a Generalized Linear model (glm) or robust linear model (rlm) to account for that.

$\endgroup$
  • $\begingroup$ Doesn't this example create negative y values which are problematic in taking the log? Also the weights of the error are much different and the lower x and y values dominate the fit because they relate to larger errors/residuals on the log scale. $\endgroup$ – Sextus Empiricus Mar 7 at 20:30
  • $\begingroup$ Right, I can get away with this since the normal error on x is quite small. Had the error been larger such that there are some negative x there are a couple different strategies to combat this although they are more based on the use case. In terms of the weights, if the normal approximation of log normal isnt sufficient, a glm model with a log normal (or gamma) error should take into account the weights (and additional dispersion) from the low values on a log scale. $\endgroup$ – badbayesian Mar 7 at 21:01
  • $\begingroup$ Instead of changing the error distribution, why not just use a power law link function (which however is not present in standard packages) or better, just perform one of the many non linear fitting methods. $\endgroup$ – Sextus Empiricus Mar 7 at 22:07
  • $\begingroup$ I have on average almost once 8000 samples a value below zero set.seed(1); y = rbeta(1e7, shape1=2,shape2=2)^(1/2.42) + rnorm(1e7, sd=.1); 1e7/sum(y<=0) has 8354.219 as output. More than half the time the lm function gives a warning message. $\endgroup$ – Sextus Empiricus Mar 7 at 22:13
  • $\begingroup$ If its only single value, you wont lose much information just dropping that value. The implementation of glm (as shown in Benjamin's post) might do a bit more than just take the logs to get around this (such as log(x+1)) but I dont know. I believe lm is giving you warnings before dropping the value. Similarly, you could try linearlize f(x) through a taylor expansion, box-cox or some other linear approximation and then running a regression. This will depend on your function and your error. This approximation method bread and butter for many econometricians. $\endgroup$ – badbayesian Mar 7 at 22:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.