Best Fit for Exponential Data

Question

I'm trying to better understand some of the theory behind fitting models that have a nonlinear link between the response and the predictors.

set.seed(1)

#Create a random exponential sample
n<-1000
y<-rexp(n=n,rate=.01)
y<-y[order(y)]
x<-seq(n)
df1<-data.frame(cbind(x,y))
plot(x,y)

Now I will attempt 4 different models fits and explain what I expect as a result versus the actual result.

m1<-lm(I(log(y))~x,data=df1)
summary(m1)
out1<-exp(predict(m1,type='response'))
lines(out1)

I expect the m1 model to be the worst fit. This is OLS with the log of the response taken before fitting the model. By taking the log of the response I have a model that, when exponentiated, cannot be negative. Therefore, the assumptions of normally distributed residuals with a constant variance cannot hold. The graph appears to under weight the tails of the distribution significantly.

Next, I will fit the model using a GLM.

m2<-glm(y~x,data=df1,family=gaussian(link='log'))
summary(m2)
out2<-predict(m2,type='response')
lines(out2,col='red')

I expect the m2 model to be of a slightly different fit than m1. This is due to the fact we are now modeling log(y+ϵ)=Xβ rather than modeling log(y)=Xβ+ϵ. I cannot justify if this fit should be better or worse from a theoretical standpoint. R gives an R-squared figure for lm() functions calls but not glm() (rather it gives and AIC score). Looking at the plot it appears that m2 is more poorly matched in the tails of the distribution than m1. The issues with normally distributed residuals should still be the same.

For m3, I relax the residual distribution assumptions by changing to a Gamma family residual distribution. Since an exponential distribution is a Gamma family distribution, I expect this model to very closely match the sampled data.

m3<-glm(y~x,data=df1,family=Gamma(link='log'))
summary(m3)
out3<-predict(m3,type='response')
lines(out3,col='blue')

AIC suggests that m3 is a better fit than m2, but it still does not appear to be a very good fit. In fact, it appears to be a worse estimate in the tails of the distribution than m1.

The m4 model will use nls just for a different approach.

m4<-nls(y~exp(a+b*x),data=df1, start = list(a = 0, b = 0))
summary(m4)
out4<-predict(m4,type='response')
lines(out4,col='yellow')

#for direct comparison
t1<-cbind(out1,out2,out3,out4,y)

This appears to very closely match m2. I also expected this to very closely fit the data.

My questions are - how come does m3 not more closely match my data? Is there a way using glm() to more closely fit this model? I realize that if I am randomly sampling then the tails of my sample data will not fit my model exactly, but these appear to be nowhere close. Perhaps letting n approach infinity would all one (or more) of the models to converge, but rerunning the code with n=100,000 actually appears to be worse fitting models (perhaps that's because with 100,000 random samples more outliers were selected and due to how the graph is present an undue amount of focus is given to these outliers?).

I also realize that an exponential distribution is not the same as exponentiating a normal distribution. That is to say, I realize that just because I took the log the of response doesn't mean that I should get a "perfect model" as a result; however, in m3 when I am fitting the models to a Gamma family distribution, I expected to be able to get a very close fit.

Thanks in advance.

Answer 1

In your simulation you are trying to model the $i$-the order statistic, $Y_{(i)}$ of an exponential distribution of an exponentially distributed sample of size $n$ as a function of $i$. This is not an exponential function, so your fitted exponential curves failed. In fact, you are essentially estimating the quantile function of the exponential distribution.

Try

curve(qexp(x/n, rate=0.01), add=TRUE)

to match your scatterplot. As a homework, derive this function and you will see that it is actually logarithmic.

If instead you want to experiment with exponentially distributed data with one predictor, use something like

y2 <- rexp(n=n, rate=exp(-1 - 0.001*x))

In this case curves 2 - 4 will pretty much overlap and reproduce the correct parameters, while curve 1 will underestimate (as expected).

Answer 2

To me it looks like your models are all off. For some reason, you decided to use the vector x = seq(n) = (1, 2, ..., n) as a covariate which doesn't make sense. Use a vector of 1s instead (an intercept-only model) and see how it goes. In fact, you don't need to create it because you can use something like:

model = glm(outcome~1, ....);

to fit an intercept-only model.

Answer 3

This is due to the fact we are now modeling log(y+ϵ)=Xβ rather than modeling log(y)=Xβ+ϵ.

My understanding is that with log link we are no longer transforming y, we are modeling yhat = e^eta, where eta is the structural model, such as a regression line. In my experience with a gnarly data set (n = 10 means, unequal sample sizes, sterr increases with mean, except for one known outlier), the StErrSq weighted gamma error model with log link gave the "best fit" by AIC, the Nweighted model gave the strongest evidence (largest LRatio). I would point out, however, that Gamma error together with StErrSq is overkill, it amounts to 2 runs at weighting by the variance, reducing the total deviance to dust.