[R] Rsquared in summary(lm)
Thomas Lumley
tlumley at u.washington.edu
Fri May 10 17:21:54 CEST 2002
On Fri, 10 May 2002, Wataru Shito wrote:
> Hi, Wouter
>
> Actually, I have a similar problem too with a simple regression.
> In my case, not only the R-square but also the estimates of intercept
> and coefficient by lm() seem different from the calculation with the
> well known formula for a simple regression.
>
> What I used is the following code. (I have just started to use R last
> week so don't blame my inmature code, please!)
> Simply,
> > ols1( y, x )
> will give you the result of the simple regression.
Can you give an example of a data set where this doesn't agree with lm()?
The code below doesn't work because the show method refers to an
apparently non-existent t.values object, but looking directly at the slots
of the value of ols1 seems to give the same results as lm
> x<-1:20
> y<-x+rnorm(20)
> summary(lm(y~x))
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-1.84289 -0.81606 0.04298 0.96198 1.57927
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.53721 0.53134 1.011 0.325
x 0.99449 0.04436 22.421 1.33e-14 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 1.144 on 18 degrees of freedom
Multiple R-Squared: 0.9654, Adjusted R-squared: 0.9635
F-statistic: 502.7 on 1 and 18 DF, p-value: 1.328e-14
> ols1(y,x)->a
> a at coefficients
$intercept
[1] 0.5372087
$coefficient
[1] 0.994491
> a at standard.errors
$intercept
[1] 0.5313436
$coeficient
[1] 0.04435571
> a at r.square
[1] 0.9654307
It would be surprising if linear regression didn't work in R in general ,
but there may be some system-specific issues
-thomas
> Wouter, could you try the following code on your data and see whether
> that's what you expect or not?
>
> I will appreciate if anyone can give me some advice why this
> differenct happens.
>
> Thankk you.
>
> Wataru Shito
>
>
> -----------------------------------------
>
> # Single Explanatory Variable Least Square Regression
> #
> library(methods)
> # create ols class
> setClass("ols", representation
> ( coefficients="list", standard.errors="list",
> r.square="numeric" ))
> setMethod("show", "ols",
> function(object)
> {
> # create row names for data.frame
> rownames <- c("(Intercept)", "X")
> # create data.frame
> z <- data.frame( row.names=rownames,
> Estimate=object at coefficients, Std.Error=object at standard.errors,
> t.value=t.values )
> cat("\n")
> print(z)
> cat( "\nR-Square:", object at r.square, "\n\n" )
> }
> )
>
> ols1 <- function( y, x ){
> size <- length(x) # number of ovservations
> xbar <- mean(x)
> ybar <- mean(y)
> Sxx <- sum( (x-xbar)^2 )
> b <- sum( (x-xbar)*(y-ybar) )/Sxx # coefficient
> a <- ybar - b*xbar # interception
> e <- y - a - b*x # residuals
> # SSE (error sum of squares)
> SSE <- sum( e^2 )
> # SST (total sum of squares)
> SST <- sum( (y-ybar)^2 )
> # SSR (regression sum of squares)
> SSR <- b^2 * Sxx
> # Coefficient of determination
> r2 <- SSR / SST
>
> # unbiased estimator of sigma^2
> s.square <- sum(e^2)/(size - 2)
> # standard error for b
> std.error.b <- sqrt( s.square/Sxx )
> # standard error for intercept
> std.error.a <- sqrt( s.square*(1/size + xbar^2/Sxx) )
> standard.errors <- list( intercept=std.error.a, coeficient=std.error.b )
> coefficients <- list( intercept=a, coefficient=b )
> new("ols", coefficients=coefficients, standard.errors=standard.errors, r.square=r2 )
> }
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>
Thomas Lumley Asst. Professor, Biostatistics
tlumley at u.washington.edu University of Washington, Seattle
^^^^^^^^^^^^^^^^^^^^^^^^
Note NEW EMAIL ADDRESS
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list