tests/OLS.r

library(regress)

## Example of ordinary least squares

x <- rep(1,15)
y <- x + rnorm(15)

mod1 <- lm(y~1)
mod2 <- regress(y~1)

## Test if coefficients are the same
if(abs(coef(mod1) - mod2$beta) > 1e-10) stop("Error with Ordinary Least Squares with constant")

## Test estimate of sigma-squared
if(abs(anova(mod1)[,3] - mod2$sigma) > 1e-10) stop("Error with Ordinary Least Squares with constant (variance)")


## Examples from lm

## Annette Dobson (1990) "An Introduction to Generalized Linear Models".
## Page 9: Plant Weight Data.  
ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2,10,20, labels=c("Ctl","Trt"))
weight <- c(ctl, trt)
lm.D9 <- lm(weight ~ group)
lm.D90 <- lm(weight ~ group - 1) # omitting intercept

rg.D9 <- regress(weight ~ group)
rg.D90 <- regress(weight ~ group - 1)

if( sum((coef(lm.D9) - rg.D9$beta)^2 + (coef(lm.D90) - rg.D90$beta)^2 > 1e-10)) stop("Error with Ordinary Least Squares")

if(abs(anova(lm.D9)[2,3] - rg.D9$sigma) > 1e-10) stop("Error with Ordinary Least Squares (variance)")
if(abs(anova(lm.D90)[2,3] - rg.D90$sigma) > 1e-10) stop("Error with Ordinary Least Squares (variance)")

Try the regress package in your browser

Any scripts or data that you put into this service are public.

regress documentation built on July 8, 2020, 6:49 p.m.