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. 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]]
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. > > 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.See ?vcov . BTW, help.search("covariance matrix") finds it.> > 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 > > >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
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. > >