[R] testing coeficients of glm
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]] __ R-help@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.
Re: [R] testing coeficients of glm
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, [EMAIL PROTECTED] 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, UKFax: +44 1865 272595 __ R-help@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.
Re: [R] testing coeficients of glm
thank you, peter On 1/24/08, Prof Brian Ripley [EMAIL PROTECTED] wrote: 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, [EMAIL PROTECTED] 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, UKFax: +44 1865 272595 -- Peter Salzman, PhD Department of Biostatistics and Computational Biology University of Rochester [[alternative HTML version deleted]] __ R-help@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.
Re: [R] testing coeficients of glm
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@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. __ R-help@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.
