[R] Bug on Mac version of lm()?
Jan de Leeuw
deleeuw at stat.ucla.edu
Sat May 11 06:57:36 CEST 2002
This is what I get with R-1.5.0 (from CRAN)
===============================================================
> source("ols1.R")
-------------------------------------
Simple regression manually calculated
-------------------------------------
(Intercept) X
Estimate -67.58065 0.97926692
Std.Error 27.91071 0.03160707
R-Square: 0.9917348
----------------------------
Simple regression using lm()
----------------------------
Call:
lm(formula = Personal.Consumption ~ Disposable.Income)
Coefficients:
(Intercept) Disposable.Income
-67.5807 0.9793
===============================================================
On Friday, May 10, 2002, at 08:37 PM, Wataru Shito wrote:
> Dear Mac users,
>
> Hi, as you might have probably read the thread of
> "[R] Rsquared in summary(lm)" on May 10, it seems that Mac version of
> lm() seem to be working incorrectly.
>
> I enclose the script to produce the result both for lm() and manual
> calculation for a simple regression. Could you run the script and
> report with the version of R, so I don't have to go through every
> builds and versions of Mac R?
>
> Here is my result where "ols1.R" is the file name for the enclosed
> script.
>
>
>> source("ols1.R")
>
> -------------------------------------
> Simple regression manually calculated
> -------------------------------------
> (Intercept) X
> Estimate -67.58065 0.97926692
> Std.Error 27.91071 0.03160707
>
> R-Square: 0.9917348
>
> ----------------------------
> Simple regression using lm()
> ----------------------------
> Call:
> lm(formula = Personal.Consumption ~ Disposable.Income)
>
> Coefficients:
> (Intercept) Disposable.Income
> 6.2334 0.9006
>
>
> This was done with Apple PowerBook G4 with OSX 10.1.4
> with XFree86 version 4.1.99.1.
>
> I got this R from CRAN archive in the binary section.
>
>> version
> _
> platform powerpc-apple-darwin5.2
> arch powerpc
> os darwin5.2
> system powerpc, darwin5.2
> status
> major 1
> minor 4.0
> year 2001
> month 12
> day 19
> language R
>
>
> Thank you for your help, in advance.
>
> ------------
> Wataru Shito
>
>
> -------------------------------------
> The contents of "ols1.R" file.
> --------------------------------------
> # Single 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")
>
> num.of.variables <- length( object at coefficients )
> # create vectors
> coeff <- array( dim=num.of.variables )
> std.err <- array( dim=num.of.variables )
> for( i in 1:num.of.variables ){
> coeff[i] <- object at coefficients[[i]]
> std.err[i] <- object at standard.errors[[i]]
> }
> # create data.frame
> z <- data.frame( row.names=rownames,
> Estimate=as.vector(coeff),
> Std.Error=as.vector(std.err) )
> cat("\n")
> print(t(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,
> coefficient=std.error.b )
> coefficients <- list( intercept=a, coefficient=b )
> new("ols", coefficients=coefficients,
> standard.errors=standard.errors, r.square=r2 )
> }
>
>
> # Demo
> #
> # "Econometric Analysis", William Greene
> # 3rd Edition, p.142, TABLE 5.1
> #
> # The result of a simple regression from the textbook is:
> # Intercept=-67.5806, b=0.979267
>
> table5.1 <- data.frame(Year = c(1970, 1971, 1972, 1973, 1974, 1975,
> 1976,
> 1977, 1978, 1979),
> Disposable.Income = c(751.6, 779.2, 810.3,
> 864.7, 857.5,
> 874.9, 906.8, 942.9, 988.8, 1015.7),
> Personal.Consumption = c(672.1, 696.8, 737.1,
> 767.9,
> 762.8, 779.4, 823.1, 864.3, 903.2, 927.6),
> row.names = c("1", "2", "3", "4", "5", "6", "7",
> "8", "9", "10"))
>
> attach( table5.1 )
> demo.ols1 <- ols1( Personal.Consumption, Disposable.Income )
> demo.lm <- lm( Personal.Consumption ~ Disposable.Income )
>
> cat( "\n-------------------------------------\n" )
> cat( "Simple regression manually calculated\n" )
> cat( "-------------------------------------" )
> print( demo.ols1 )
>
> cat( "----------------------------\n" )
> cat( "Simple regression using lm()\n" )
> cat( "----------------------------" )
> print( demo.lm )
> detach( table5.1 )
>
>
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> .-.-.-
> 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
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
> ._._._
>
>
===
Jan de Leeuw; Professor and Chair, UCLA Department of Statistics;
US mail: 9432 Boelter Hall, Box 951554, Los Angeles, CA 90095-1554
phone (310)-825-9550; fax (310)-206-5658; email: deleeuw at stat.ucla.edu
homepage: http://www.stat.ucla.edu/~deleeuw
========================================================
No matter where you go, there you are. --- Buckaroo Banzai
http://www.stat.ucla.edu/~deleeuw/sounds/nomatter.au
========================================================
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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