[R] Test if beta is different from something other than 0

Thu Jan 10 18:31:58 CET 2002

 >  -----Original Message-----
 >  From: Moffet, Corey [mailto:Corey.Moffet at ttu.edu]
 >  Sent: Thursday, January 10, 2002 8:47 AM
 >  To: 'r-help at stat.math.ethz.ch'
 >  Subject: [R] Test if beta is different from something other than 0
 >  
 >  
 >  Is there a function/package that will allow you to test the 
 >  hypothesis beta1
 >  = x in a
 >  simple linear regression, where x is a constant?  

I've written a funcgion 'glh.test' for testing general linear hypotheses.  

For this case, you can use glh.test as:
	> x <- rnorm(100)
	> y <- rnorm(100,x+1)
	> reg <- lm(y~x)
	> # test the hypothesis that (b1*0, b1*1) == 1
	> glh.test( reg,            c(   0,    1), c(1) )

		 Test of General Linear Hypothesis 	
	Call:
	glh.test(reg = reg, cm = c(0, 1), d = c(1))
	F = 2.331, df1 =  1, df2 = 98, p-value = 0.1300 

I'm pasting the function and the documentation below. 'glh.test' will be
part of the next release of the gregmisc library, which will be out as soon
as I can get it through the publication review process.

-Greg

##################### start of glh.test.R ############################
# $Id: glh.test.R,v 1.3 2001/12/19 20:05:27 warneg Exp $
#
# $Log: glh.test.R,v $
# Revision 1.3  2001/12/19 20:05:27  warneg
#
# - Removed extra element of return object.
#
# Revision 1.2  2001/12/18 21:34:25  warneg
# - Modified to work correctly when obj is of class 'aov' by specifying
#   summary.lm instead of summary.  This ensures that the summary object
#   has the fields we need.
#
# - Moved detailed reporting of results from 'print' to 'summary'
#   function and added a simpler report to 'print'
#
# Revision 1.1  2001/12/18 00:43:23  warneg
#
# Initial checkin.
#
#

glh.test <- function( reg, cm, d=rep(0, nrow(cm)) )
{

  if( !is.matrix(cm) && !is.data.frame(cm) )
    cm <- matrix(cm, nrow=1)

  if ( !( "lm" %in% class(reg) ) )
    stop("Only defined for lm,glm objects")

  bhat <- summary.lm(reg)$coefficients[,1,drop=F]
  XpX <- summary.lm(reg)$cov.unscaled
  df <- reg$df.residual
  msr <- summary.lm(reg)$sigma  # == SSE / (n-p)
  r <- nrow(cm)

  if ( ncol(cm) != length(bhat) ) stop(  
                   paste( "\n Dimension of ",
                         deparse( substitute( cm ) ), ": ",
o                         paste( dim(cm), collapse="x" ),
                         ", not compatible with no of parameters in ",
                         deparse( substitute( reg ) ), ": ",
                         length(bhat), sep="" ) )

  #                        -1
  #     (Cb - d)' ( C (X'X)   C' ) (Cb - d) / r
  # F = ---------------------------------------
  #                 SSE / (n-p)
  #

  Fstat <- t(cm %*% bhat - d) %*% solve((cm %*% XpX %*% t(cm))) %*% (cm %*%
bhat - d) / r / msr^2 

  p <- 1-pf(Fstat,r,df)

  retval <- list()
  retval$call <- match.call()
  retval$statistic <- c(F=Fstat)
  retval$parameter <- c(df1=r,df2=df)
  retval$p.value <- p
  retval$conf.int <- NULL
  retval$estimate <- cm%*%bhat
  retval$null.value <- d
  retval$method <- "Test of General Linear Hypothesis"
  retval$data.name <- deparse(substitute(reg))
  retval$matrix <- cm
  colnames(retval$matrix) <- names(reg$coef)

  class(retval) <- c("glh.test","htest")

  retval
}

print.glh.test <- function(x, digits = 4 )
{
    cat("\n")
    cat("\t",x$method, prefix = "\t")
    cat("\n")
    cat("Call:\n")
    print(x$call)

    if (!is.null(x$statistic)) 
        cat(names(x$statistic), " = ", format(round(x$statistic, 
            4)), ", ", sep = "")
    if (!is.null(x$parameter)) 
        cat(paste(names(x$parameter), " = ", format(round(x$parameter, 
            3)), ",", sep = ""), "")
    cat("p-value =",
        format.pval(x$p.value, digits = digits), 
        "\n")
    cat("\n")
  }

summary.glh.test <- function(x, digits = 4 )
{
    cat("\n")
    cat("\t",x$method, prefix = "\t")
    cat("\n")
    cat("Regression: ", x$data.name, "\n")
    cat("\n")
    cat("Null Hypothesis: C %*% Beta-hat = d \n")
    cat("\n")
    cat("C matrix: \n")
    print(x$matrix, digits=digits)
    cat("\n")
    cat("d vector: \n")
    print(x$null.value, digits=digits)
    cat("\n")
    cat("C %*% Beta-hat: \n")
    print(c(x$estimate))
    cat("\n")

    if (!is.null(x$statistic)) 
        cat(names(x$statistic), " = ", format(round(x$statistic, 
            4)), ", ", sep = "")
    if (!is.null(x$parameter)) 
        cat(paste(names(x$parameter), " = ", format(round(x$parameter, 
            3)), ",", sep = ""), "")
    cat("p-value =",
        format.pval(x$p.value, digits = digits), 
        "\n")
    cat("\n")
  }
##################### end of glh.test.R ###########################

##################### start of glh.test documentation
###########################
glh.test              package:gregmisc              R Documentation
Test a General Linear Hypothesis for a Regression Model

Description:

     Test, print, or summarize a general linear hypothesis for a
     regression model

Usage:

     glh.test(reg, cm, d=rep(0, nrow(cm)) )
     print.glh.test(x, digits=4)
     summary.glh.test(x, digits=4)

Arguments:

     reg: Regression model 

      cm: matrix .  Each row specifies a linear combination of the
          coefficients 

       d: vector specifying the null hypothis values for each linear
          combination

       x: glh.test object

  digits: number of digits.

Details:

     Test the general lineary hypothesis C %*% beta_hat == d  for the
     regression model `reg'.

     The test statistic is obtained from the formula: 

 F = (C Beta-hat - d)' ( C (X'X)^-1   C' ) (C Beta-hat - d) / r / ( SSE /
(n-p) )

     Under the null hypothesis, f will follow a F-distribution with r
     and n-p degrees of freedom.

Value:

     Object of class `c("glh.test","htest")' with elements: 

   call : Function call that created the object

statistic : F statistic

parameter: vector containing the numerator (r) and denominator (n-p)
          degrees of freedom 

 p.value: p-value

estimate: computed estimate for each row of `cm' 

null.value: d 

  method: description of the method 

data.name: name of the model given for `reg' 

  matrix: matrix specifying the general linear hypotheis (`cm')

Author(s):

     Gregory R. Warnes gregory_r_warnes\@groton.pfizer.com

References:

     R.H. Myers, Classical and Modern Regression with Applications, 
     2nd Ed, 1990, p. 105