[R] Exponential regression (Y = exp(a*X)) and standard error of Ŷi

David Winsemius dwinsemius at comcast.net
Tue Oct 28 02:44:00 CET 2008


I am reasonably sure that the log link with a zero intercept  
duplicates your efforts with nls().
Take a look at this:

 > txt <- "X         Y
+ 2         2.927064
+ 5         14.6582
+ 4         8.567706
+ 3         5.007817
+ 1         1.710867
+ 6         25.07823
+ 4         8.567706
+ 7         42.9055
+ 2         2.927064
+ 8         125.5872
+ 8         125.5872
+ 7         42.9055    "

 > dta <- read.table(textConnection(txt),header=TRUE)

 > log.lin <- glm(Y~X, data=dta, family=quasi(link="log",  
variance="constant"))
 > log.lin

Call:  glm(formula = Y ~ X, family = quasi(link = "log", variance =  
"constant"),      data = dta)

Coefficients:
(Intercept)            X
     -2.6731       0.9368

Degrees of Freedom: 11 Total (i.e. Null);  10 Residual
Null Deviance:	    22500
Residual Deviance: 251.1 	AIC: NA

### so not terribly close to what you were seeing
### You never did say what you got from R

 > X <- dta$X
 > Y <- dta$Y
 > exp1 <- nls(Y~exp(X*a), start=list(a=0.3), trace=T)
29193.26 :  0.3
2973.946 :  0.5594376
1598.388 :  0.577918
1209.806 :  0.5929772
1208.918 :  0.592337
1208.918 :  0.5923281
 > log.lin <- glm(Y~X+0, data=dta, family=quasi(link="log",  
variance="constant"))
 > log.lin

Call:  glm(formula = Y ~ X + 0, family = quasi(link = "log", variance  
= "constant"),      data = dta)

Coefficients:
      X
0.5923

### So a log link and constant variance model produces a parameter  
estimate that matches the nls call:

Degrees of Freedom: 12 Total (i.e. Null);  11 Residual
Null Deviance:	    35460
Residual Deviance: 1209 	AIC: NA

 > summary(log.lin)

Call:
glm(formula = Y ~ X + 0, family = quasi(link = "log", variance =  
"constant"),
     data = dta)

Deviance Residuals:
      Min        1Q    Median        3Q       Max
-20.2940   -5.9720   -1.5133   -0.2812   11.3103

Coefficients:
   Estimate Std. Error t value Pr(>|t|)
X 0.592328   0.007199   82.28   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’  
0.1 ‘ ’ 1

(Dispersion parameter for quasi family taken to be 109.9048)

     Null deviance: 35460.9  on 12  degrees of freedom
Residual deviance:  1208.9  on 11  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4
###-------------

The std error is part of the output of the summary method for glm  
objects. Perhaps if you think about your fitting process in terms of  
generalized linear models and read a bit more about what SPSS is  
really doing, you will have a better understanding of the issues  
involved (choice of link and variance model).

HTH;
David Winsemius
Heritage Labs

On Oct 27, 2008, at 11:05 AM, alemu gonsamo wrote:

> r-help at lists.R-project.org
>
> Hello
>
> First I want to implement exponential regression in R, with out  
> constant for the following formula.
>
> Y = exp(a*X)
>
> ‘a’ is coefficient I wanted to determine. That I could do also in  
> SPSS but my question is rather to estimate the ‘standard error of   
> Ŷi  at each Xi. This is called in SPSS ‘satndard error of mean  
> prediction’ or generally known for non-linear regression as  
> ‘asymptotic standard error’. This is different from residual.
>
> Below is the example data set for which I wanted to calculate the  
> coefficient ‘a’ for the exponential regression in the form stated  
> above and ‘standard error of Ŷi’. For this specific data,  
> ‘a’ is computed using SPSS and the result is 0.5620. So, Y =  
> exp(0.5620*X).
>
>
> X         Y                 SE of Ŷi
> 2         2.927064         ?
> 5         14.6582           ?
> 4         8.567706         ?
> 3         5.007817         ?
> 1         1.710867         ?
> 6         25.07823         ?
> 4         8.567706         ?
> 7         42.9055           ?
> 2         2.927064         ?
> 8         125.5872         ?
> 8         125.5872         ?
> 7         42.9055          ?
>
> My questions are (1), how I implement exponential regression in R  
> and (2) how to calculate the ‘standard error of Ŷi’’ for each  
> Xi.
>
> Well,
>> exp1 <- nls(Y~exp(X*a), start=list(a=0.3), trace=T)
>
> from this I can get summary statistics, 'a' is different from what i  
> obtained using SPSS. Can someone help please atleast with getting SE  
> of Ŷi.
>
> I wrote very long to make the question clear.
>
> thanks in advance
>
> ggalex
>
> Helsinki, Finland
> ggalex2002 at yahoo.com
>
>
>
>
>
> ______________________________________________
> 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