sigr: Simple Significance Reporting

[This article was first published on R – Win-Vector Blog, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

sigr is a simple R package that conveniently formats a few statistics and their significance tests. This allows the analyst to use the correct test no matter what modeling package or procedure they use.

Sigr

Model Example

Let’s take as our example the following linear relation between x and y:

library('sigr')
set.seed(353525)
d <- data.frame(x= rnorm(5))
d$y <- 2*d$x + 0.5*rnorm(nrow(d))

stats::lm() has among the most complete summaries of all models in R, so we easily get can see the quality of fit and its significance:

model <- lm(y~x, d=d)
summary(model)
## 
## Call:
## lm(formula = y ~ x, data = d)
## 
## Residuals:
##       1       2       3       4       5 
##  0.3292  0.3421  0.0344 -0.1671 -0.5387 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   0.4201     0.2933   1.432  0.24747   
## x             2.1117     0.2996   7.048  0.00587 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4261 on 3 degrees of freedom
## Multiple R-squared:  0.943,  Adjusted R-squared:  0.9241 
## F-statistic: 49.67 on 1 and 3 DF,  p-value: 0.005871

sigr::wrapFTest() can render the relevant model quality summary.

cat(render(wrapFTest(model),
    pLargeCutoff=1))

F Test summary: (R2=0.94, F(1,3)=50, p=0.0059).

Note: sigr reports the un-adjusted R-squared, as this is the one that significance is reported for in summary::lm().

sigr also carries around the important summary components for use in code.

unclass(wrapFTest(model))
## $test
## [1] "F Test"
## 
## $numdf
## [1] 1
## 
## $dendf
## [1] 3
## 
## $FValue
## [1] 49.66886
## 
## $R2
## [1] 0.9430403
## 
## $pValue
## [1] 0.005871236

In this function sigr is much like broom::glance() or modelr::rsquare().

broom::glance(model)
##   r.squared adj.r.squared     sigma statistic     p.value df    logLik
## 1 0.9430403     0.9240538 0.4261232  49.66886 0.005871236  2 -1.552495
##       AIC      BIC deviance df.residual
## 1 9.10499 7.933304 0.544743           3
modelr::rsquare(model, d)
## [1] 0.9430403

This is something like what is reported by caret::train() (where caret controls the model training process).

cr <- caret::train(y~x, d, 
                   method = 'lm')
## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info =
## trainInfo, : There were missing values in resampled performance measures.
cr$results
##   intercept      RMSE  Rsquared    RMSESD   RsquaredSD
## 1      TRUE 0.9631662 0.9998162 0.6239989 0.0007577834
summary(cr$finalModel)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      X1      X2      X3      X4      X5 
##  0.3292  0.3421  0.0344 -0.1671 -0.5387 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   0.4201     0.2933   1.432  0.24747   
## x             2.1117     0.2996   7.048  0.00587 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4261 on 3 degrees of freedom
## Multiple R-squared:  0.943,  Adjusted R-squared:  0.9241 
## F-statistic: 49.67 on 1 and 3 DF,  p-value: 0.005871

(I presume cr$results$Rsquared is a model quality report and not a consistency of cross-validation procedure report. If it is a model quality report it is somehow inflated, perhaps by the resampling procedure. So I apologize for using either caret::train() or its results incorrectly.)

Data example

The issues in taking summary statistics (and significances) from models include:

  • Working only from models limits you to models that include the statistic you want.
  • Working only from models is mostly "in-sample." You need additional procedures for test or hold-out data.

With sigr it is also easy to reconstruct quality and significance from the predictions, no matter where they came from (without needing the model data structures).

In-sample example

d$pred <- predict(model, newdata = d)
cat(render(wrapFTest(d, 'pred', 'y'),
    pLargeCutoff=1))

F Test summary: (R2=0.94, F(1,3)=50, p=0.0059).

Notice we reconstruct the summary statistic and significance, independent of the model data structures. This means the test is generic and can be used on any regression (modulo informing the significance model of the appropriate number of parameters). It also can be used on held-out or test data.

In this mode it is a lot like ModelMetrics::rmse() or caret::postResample().

ModelMetrics::rmse(d$y, d$pred)
## [1] 0.3300736
caret::postResample(d$y, d$pred)
##      RMSE  Rsquared 
## 0.3300736 0.9430403

Out of sample example

If we had more data we can look at the quality of the prediction on that data (though you have to take care in understanding the number of degrees of freedom is different for held-out data).

d2 <- data.frame(x= rnorm(5))
d2$y <- 2*d2$x + 0.5*rnorm(nrow(d2))
d2$pred <-  predict(model, newdata = d2)
cat(render(wrapFTest(d2, 'pred', 'y'),
    pLargeCutoff=1))

F Test summary: (R2=0.76, F(1,3)=9.6, p=0.054).

Exotic model example

We could have used glmnet instead of lm:

library("glmnet")
d$one <- 1 # make sure we have at least 2 columns
xmat <- as.matrix(d[, c('one','x')])
glmnetModel <- cv.glmnet(xmat, d$y)
## Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations
## per fold
summary(glmnetModel)
##            Length Class  Mode     
## lambda     53     -none- numeric  
## cvm        53     -none- numeric  
## cvsd       53     -none- numeric  
## cvup       53     -none- numeric  
## cvlo       53     -none- numeric  
## nzero      53     -none- numeric  
## name        1     -none- character
## glmnet.fit 12     elnet  list     
## lambda.min  1     -none- numeric  
## lambda.1se  1     -none- numeric
d$predg <- as.numeric(predict(glmnetModel, 
                              newx= xmat))
cat(render(wrapFTest(d, 'predg', 'y'),
    pLargeCutoff=1))

F Test summary: (R2=0.91, F(1,3)=30, p=0.012).

Plotting

Because sigr can render to "LaTex" it can (when used in conjunction with latex2exp) also produce formatted titles for plots.

library("ggplot2")
library("latex2exp")


f <- paste0(format(model$coefficients['x'], digits= 3), 
            '*x + ',
            format(model$coefficients['(Intercept)'], digits= 3))
title <- paste0("linear y ~ ", f, " relation")
subtitle <- latex2exp::TeX(render(wrapFTest(d, 'pred', 'y'), 
                                          format= 'latex'))
ggplot(data=d, mapping=aes(x=pred, y=y)) + 
  geom_point() + geom_abline(color='blue') +
  xlab(f) +
  ggtitle(title, 
          subtitle= subtitle)
NewImage

Conclusion

sigr computes a few statistics or summaries (with standard significance estimates) and returns the calculation in both machine readable and formatted forms. For non-standard summaries sigr includes some resampling methods for significance estimation. For formatting sigr tries to get near the formats specified by "The American Psychological Association." sigr works with models (such as lm and glm(family='binomial')) or data, and is a small package with few dependencies.

To leave a comment for the author, please follow the link and comment on their blog: R – Win-Vector Blog.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)