6

I come from an SPSS background and am attempting to move to R for it's superior flexibility and data manipulation abilities. I have some concerns however as to whether the lm() is really using partial correlations.

I'm basically trying to run a linear regression, using something similar to the "enter" setting in SPSS, which essentially builds the model one variable at a time, reporting the change in $R^2$ with each additional variable. This allows you to determine how much predictive power each variable adds to the model.

When I run the same analysis in R however, I don't get any information on the $R^2$ contributed by individual variables, and I'm not even sure that it's using partial corrrelations to calculate the p-values it's reporting!

My code follows:

summary(m1 <- lm(totalprop ~ cos(Angle) + Alignment + colour + 
  Angle*Alignment, dataset))

My questions:

  1. Does R use partial correlations to determine reported p-values from lm()?
  2. How can I make R report change in $R^2$ with each additional variable.
  3. How can I make R act like SPSS in calulating the model piece by piece? Is this possible without running multiple iterations of the lm() function? If not, how does one control for the effects of covariates in R?
analystic
  • 725
  • 1
  • 5
  • 21
  • I've just added the final part of my question, as I don't think it was clear - I want R to be SPSS for me for a moment! I'm aware it's bad practise to try to fit a foot into a glove this way, but I assume there must be a way to do standard statistical control without running iterations of models on residuals as described here: http://stats.stackexchange.com/questions/17336/how-exactly-does-one-control-for-other-variables – analystic Sep 22 '12 at 12:37

3 Answers3

6

Another function that you might find useful is lm.sumSquares in the lmSupport package. Basically, if you have the following model:

mod1 <- lm(dv ~ iv1 + iv2 + iv3, data = d)

Executing the following will give you the delta R squared (the semipartial correlation squared) and PRE (the partial correlation squared) for iv1, iv2, and iv3:

lm.sumSquares(mod1)

lmSupport also contains the function lm.deltaR2 that allows you to compare custom models against each other to obtain the change in R squared and the F-statistic associated with the R squared change. So, going back to the above example, if you also had this second model:

mod2 <- lm(dv ~ iv1 + iv2 + iv3 + iv4 + iv5, data = d)

Then you could do the following to obtain the R squared change for adding iv4 and iv5 and the accompanying F statistic and p-value:

lm.deltaR2(mod1, mod2)

Hope that helps!

Patrick S. Forscher
  • 3,122
  • 23
  • 43
  • I'm marking this as the answer, although I haven't yet run the code. When I get a chance I may revise that, but based on your description this looks like exactly what I need. Hopefully, it'll allow me to run several models and watch how the partial correlations change as I add variables, I can then fudge a hierarchical regression from that with some luck! – analystic Sep 26 '12 at 13:59
  • Finally had a chance to have a look at this. Looks almost perfect for my needs, but doesn't include the p values for those eta2 values. I'm going to leave it as the answer as I didn't originally ask for a package that gave the p values on the partial R2! I'm guessing maybe when my R skills get good enough I can modify the package to display this info as well, as it already includes partial SS and SSR. – analystic Sep 28 '12 at 14:59
  • You can find the p-values for the eta squared / delta R squared values in the plain summary output for lm. If you're testing a set of variables, lm.deltaR2 gives p-values for the delta R squared it gives you. – Patrick S. Forscher Sep 28 '12 at 15:20
  • Hi Patrick, I think the p values given in the regular lm output are for the zero order correlation, not for the partial correlation? – analystic Sep 29 '12 at 11:46
  • Sorry, I see what you mean now, but lm.deltaR2() gives the overall significance of the difference in explanatory power between two models, but what I'm looking for is the significance of each variable after the variance of the other variables has been explained. I think I need to be clearer on the difference between overall R2 explained by the model and the R2 explained by each individual variable, after the effect of all other variables has already been partialled out... – analystic Sep 29 '12 at 11:53
  • 1
    For those reading this currently, lm.deltaR2 has been deprecated, use modelCompare instead. – Rilcon42 Mar 17 '16 at 15:36
5

R doesn't use partial correlations to determine reported p-values from lm.

Indeed, you can examine the code that computes the p-values for yourself; if you type 

> summary.lm

you will see the code of the function includes the following:

    ans$coefficients <- cbind(est, se, tval, 2 * pt(abs(tval), 
    rdf, lower.tail = FALSE))

which shows that for each variable, the p-value is computed from the t value by looking up t tail-areas

It is quite easy to compute changes in R^2 between models if you need to do that, but I expect one of the many packages on CRAN does it already.

Glen_b
  • 257,508
  • 32
  • 553
  • 939
  • Thanks for the info Glen; is there any chance you could be more specific about how I can get the answers I'm looking for from R? Is there a specific package that allows evaluation of the change in R^2 without iterating through several models like rolando2 suggests above? This approach seems to violate the DRY principle pretty badly. Alternatively, is there a package that allows me to run a multiple regression using partial correlations rather than the zero order correlations, and assess the significance of those partial correlations? – analystic Sep 22 '12 at 14:42
5

I'm no R expert but I've looked into this same issue (your 2nd and 3rd questions) and found the following.

  1. You do need first to run the lm command for each model separately.
  2. If a formal test of the F associated with the change is of interest to you, you can formally compare successive models using anova, e.g., anova(MyModel1, MyModel2).
  3. This anova test will only work if each solution is based on the same subset of cases. So you may need to run a command at the outset such as mydata=na.omit(mydata).

EDIT: the relaimpo package's calc.relimp command ("relative importance") calculates delta rsq for each predictor when entered last (its semipartial rsq) and/or when entered first (its zero-order rsq).

install.packages("relaimpo", dep=c("Depends"))  
library(relaimpo)  
calc.relimp( MyRegressionModel, type = c("last", "first") )
rolando2
  • 11,645
  • 1
  • 39
  • 60
  • This is a helpful answer and a suggestion I've come across previously... As the link I posted above suggested, this would also allow you to calculate a partial correlation if you regressed each iteration of the model on the residuals of the model before it. However I still think it violates the DRY principle a little, as I'd still have to build a succession of models, and still wouldn't be able to view the zero order correlations against the partials... – analystic Sep 22 '12 at 15:13
  • @analystic What is the DRY principle? – mark999 Sep 22 '12 at 21:03
  • Don't Repeat Yourself :) If I have to build the model one step at a time, that means I have to repeatedly call the same function (`lm`), and calculate the difference in R^2 between iterations of the model. – analystic Sep 23 '12 at 03:16
  • 1
    @analystic It's relatively easy to avoid repeating yourself in R - you can, for example, loop through model updates (by calling update) and keep results in a list. There are even fancier ways to do things. Any time you think you're repeating yourself, your probably failing to use (or write) some function or vectorization that would take care of it (the *apply functions, for example, often avoid repeating something) – Glen_b Sep 27 '12 at 02:57
  • @Glen That's true, however, I did think there must have been a pre built package to do this task quickly, as it's a very standard part of regression analysis. To be honest, I'm surprised it's not built into the R language's lm function from the start! – analystic Sep 28 '12 at 15:02
  • 1
    @analystic I also think there's likely a pre built package to do this task. I don't know what that might be. If I wanted to do something like that I imagine I'd write a function. For myself, I'm not at all surprised it's not in the lm function. – Glen_b Sep 29 '12 at 01:22