[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