[R] Bug on Mac version of lm()?

[Wataru Shito]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._