[R] testing coeficients of glm

[Achim Zeileis]

On Thu, 24 Jan 2008, peter salzman wrote:

> Dear list,
>
> i'm trying to test if a linear combination of coefficients of glm is equal
> to 0. For example :
> class 'cl' has 3 levels (1,2,3) and 'y' is a response variable. We want to
> test  H0: mu1 + mu2 - mu3 =0 where mu1,mu2, and mu3 are  the means for each
> level.

Look at
  linear.hypothesis()
in package "car".

## data (reduce treatment to three levels)
data("OrchardSprays")
os <- subset(OrchardSprays, treatment %in% c("A", "B", "C"))
os$treatment <- factor(os$treatment)

## model and test
fm <- lm(decrease ~ 0 + treatment, data = os)
coef(fm)
linear.hypothesis(fm, "treatmentA + treatmentB - treatmentC = 0")

See the accompanying documentation for an overview which extractor
functions (such as coef, vcov, etc.) are used to compute the test.

hth,
Z

> for me, the question is how to get the covariance matrix of the estimated
> parameters from glm. but perhaps there is a direct solution in one of the
> packages.
>
> i know how to solve this particular problem (i wrote it below) but i'm
> curious about the covariance matrix of coefficient as it seems to be
> important.
>
> the R code example :
> ###
> nObs <- 10
> cl <- as.factor( sample(c(1,2,3),nObs,replace=TRUE) )
> y <- rnorm(nObs)
>
> model <- glm(y ~ cl)
> b <- model$coefficients
> H <- c(1,1,-1)  # we want to test H0: Hb = 0
>
> ### the following code will NOT run unless we can compute covModelCoeffs
>
> #the mean of Hb is
> mu = H %*% model$coefficients
> #the variance is HB is
> var = H %*% covModelCoeffs %*% t(H)
>
> p.val <- 2 * pnorm( -abs(mu), mean=0, sd=sqrt(var),lower.tail = TRUE)
>
>
> how do i get the covariance matrix of the estimated parameters ?
>
> thanks,
> peter
>
> P.S. the simple solution for this particular problem:
>
> ## get the mean for each level
>  muV <- by(y,cl,mean)
> ## get the variance for each level
>  varV <- by(y,cl,var)
>
> ## the mean of Hb is
> muHb <- H %*% muV
> ## because of independence, the variance of Hb is
> varHb <- sum(varV)
>
> ## the probability of error, so-called p-value:
> p.val <- 2 * pnorm( -abs(muHb), mean=0, sd=sqrt(varHb),lower.tail = TRUE)
>
> thanks again,
> peter
>
>
> --
> Peter Salzman, PhD
> Department of Biostatistics and Computational Biology
> University of Rochester
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>