5

I'm wondering if there's an easy way of calculating an F-statistic / p-value for a subset of model coefficients. Particularly in R? I'm not sure what test would be needed to calculate this. For example,

summary(lm(a~w+x+y+z))

will give me the F-statistic and p-value for the whole model, but is it possible to extract the F-statistic and p-value for just the terms y and z?

A faster way of doign something like this:

matrixOfResponses <- cbind(c(1,2,3,4,5), c(4,3,2,4,5), c(5,3,2,23,4), c(1,2,4,3,1,))
pValsOut <- numeric()
for(i in 1:ncol(matrixOfResponses))
{
pValsOut[i] <- anova(lm(matrixOfResponses[,i]~mMat), lm(matrixOfResponses[,i]~mMatReduced))$'Pr(>F)'[2]
}

So basically something like the above becomes very slow when "matrixOfRespones" contains a huge number of variables, for example a huge number of gene expression levels.

ScanSally
  • 71
  • 1
  • 4
  • Welcome to the site, @ScanSally. This appears to be *only* a question about how to get R to do this for you. If that's true, this Q would be off-topic for CV (see our [FAQ](http://stats.stackexchange.com/faq)), but on-topic for [Stack Overflow](http://stackoverflow.com/). Are you also interested in understanding (eg) the use of the F-test as a simultaneous test to assess multiple variables, or some other substantive statistical issue? If so, please edit to verify, if not, flag your Q & we'll migrate it for you (*please don't cross-post, though*). – gung - Reinstate Monica Feb 25 '13 at 20:35
  • It was more a general stats question, so I'll edit the question! – ScanSally Feb 25 '13 at 23:39
  • The updated portion of the question really is just a programming question. The answer is probably *vectorization*. You might want to start by looking here: [vectorize-my-thinking-vector-operations-in-r](http://stackoverflow.com/questions/445059/). – gung - Reinstate Monica Mar 06 '13 at 18:28

2 Answers2

2

I suspect what you are interested in is testing both $y$ and $z$ simultaneously. (NB, if not, @jebyrnes' answer addresses your question.) To do this, you fit two models and assess them with a nested model test (often called "$F$ change test", or "$R^2$ change test"): 

fullMod    = lm(a~w+x+y+z)
reducedMod = lm(a~w+x)
anova(fullMod, reducedMod)

For more conceptual understanding, my answer here Testing for moderation with continuous vs categorical moderators talks about the $F$ change test (albeit in a different context).

Community
  • 1
gung - Reinstate Monica
  • 132,789
  • 81
  • 357
  • 650
  • Great this is exactly what I was looking for. Thanks!! – ScanSally Feb 25 '13 at 23:40
  • That's great, @ScanSally. If you want to know more, I recommend reading my other answer, a lot of it is about testing multiple parameters when they are indicator variables used in an interaction, but the principle is the same. – gung - Reinstate Monica Feb 26 '13 at 00:02
  • This may be a long shot, but is there any chance its possible to do something like the above for a large number of response variables. I.e. "a" is a matrix of ~50,000 different response variables, so fullMod and reducedMod will be objects of the class c("mlm" "lm"). It seems really slow to do this using a for loop, but maybe there is a much faster way that I just don't know about? – ScanSally Mar 06 '13 at 17:15
  • In theory, it's certainly possible to test a nested model w/ just a few variables against a full model w/ 50k variables. I don't know anything about what performance issues that might exist there in practice though. If you are asking about the possibility of scrolling through large numbers of variables to find out which are 'the right ones', you should read my answer here: [algorithms-for-automatic-model-selection](http://stats.stackexchange.com/questions/20836//20856#20856). – gung - Reinstate Monica Mar 06 '13 at 17:24
  • Thanks for the answer, but I think I mean something a little different! I'm not trying to do variable selection, just save the trouble of having to do a for loop to compare fullMod and reducedMod for different "a" (i.e. different response variables) above. – ScanSally Mar 06 '13 at 17:39
  • Sorry, I misread your question. I gather you want to compare a model w/ multiple **response** variables to a model w/ just 1. I don't think that would be valid, but it's an interesting question. – gung - Reinstate Monica Mar 06 '13 at 17:54
  • No not really, I just want a faster way of doing something like this: Oh wait I can put code here, hang on I'll write it above in the main text. – ScanSally Mar 06 '13 at 17:57
1

That summary is in and of itself an object. So, 

tab <- summary(lm(a~w+x+y+z))
tab$coefficients

Will give you the coefficient table. It's just a matrix. You can then extract the appropriate entries from there. anova and Anova would similarly.

jebyrnes
  • 835
  • 6
  • 16