You’ll learn about these later.

## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union


My Favorite R Packages

You will want these for reference. I tried to include only those that were not part of in class assignments.

e1071

e1071 calculates skewness. The formula it uses depends on the type assigned, so refer to this page and choose the one that most reflects your class:

https://www.rdocumentation.org/packages/e1071/versions/1.7-17/topics/skewness

library(e1071)
skewness(mtcars$mpg, type = 1)
## [1] 0.6404399
skewness(mtcars$mpg, type = 2)
## [1] 0.6723771
skewness(mtcars$mpg, type = 3)
## [1] 0.610655


psych

psych has the function describe. This automatically generates a table of descriptive statistics. You can filter out what you don’t want.

library(psych)
desc <- describe(mtcars, quant=c(.25, .75))
desc$term <- rownames(desc) #without this, your variables won't be included in the nice table function, but you also can't print the raw df because the rownames are not numbers.

desc.final <- desc %>%
  select(term, min, Q0.25, mean, Q0.75, max, sd, skew)

colnames(desc.final) <- c(
  "Variable",
  "Minimum",
  "25th Percentile",
  "Mean",
  "75th Percentile",
  "Maximum",
  "Standard Deviation",
  "Skew"
)


Rempsyc

Rempsyc has the function nice_table. This automatically generates a publishable table from something like your descriptive statistics. I usually partner this with psych.

library(rempsyc)
## Suggested APA citation: Thériault, R. (2023). rempsyc: Convenience functions for psychology. 
## Journal of Open Source Software, 8(87), 5466. https://doi.org/10.21105/joss.05466
nice_table(desc.final,
           title = c("Table 1: Descriptive Statistics of Variables"))

Table 1: Descriptive Statistics of Variables

Variable

Minimum

25th Percentile

Mean

75th Percentile

Maximum

Standard Deviation

Skew

mpg

10.40

15.43

20.09

22.80

33.90

6.03

0.61

cyl

4.00

4.00

6.19

8.00

8.00

1.79

-0.17

disp

71.10

120.83

230.72

326.00

472.00

123.94

0.38

hp

52.00

96.50

146.69

180.00

335.00

68.56

0.73

drat

2.76

3.08

3.60

3.92

4.93

0.53

0.27

wt

1.51

2.58

3.22

3.61

5.42

0.98

0.42

qsec

14.50

16.89

17.85

18.90

22.90

1.79

0.37

vs

0.00

0.00

0.44

1.00

1.00

0.50

0.24

am

0.00

0.00

0.41

1.00

1.00

0.50

0.36

gear

3.00

3.00

3.69

4.00

5.00

0.74

0.53

carb

1.00

2.00

2.81

4.00

8.00

1.62

1.05


Broom

nice_table can also be used for your models, but the output when you generate a linear model is pretty messy.

car_lm <- lm(mpg ~ cyl + disp + hp + drat + wt, data = mtcars)
car_lm
## 
## Call:
## lm(formula = mpg ~ cyl + disp + hp + drat + wt, data = mtcars)
## 
## Coefficients:
## (Intercept)          cyl         disp           hp         drat           wt  
##    36.00836     -1.10749      0.01236     -0.02402      0.95221     -3.67329


To make it usable for nice table, we use the tidy() function from the broom package.

library(broom)
car_lm %>% tidy() %>% nice_table(title = "Linear Model Table")

Linear Model Table

Term

estimate

std.error

statistic

p

(Intercept)

36.01

7.57

4.76

< .001***

cyl

-1.11

0.72

-1.55

.134

disp

0.01

0.01

1.04

.308

hp

-0.02

0.01

-1.81

.082

drat

0.95

1.39

0.68

.500

wt

-3.67

1.06

-3.47

.002**


Leaps

Leaps lets us view the Mallow’s CP, BIC, and Adjusted R Squared for model selection.

library(leaps)
car_sub <- regsubsets(mpg ~ cyl + disp + hp + drat + wt, data = mtcars)
summary(car_sub)$cp
## [1] 15.221274  3.687548  3.427820  4.468710  6.000000
summary(car_sub)$bic
## [1] -37.79462 -46.34824 -45.41594 -43.08923 -40.19523
summary(car_sub)$adjr
## [1] 0.7445939 0.8185189 0.8263446 0.8262103 0.8227219


Performance

Sometimes you’ll be building a logistic regression model, and you’ll want to look at a plot of your residuals to check for fit. But oh no! Your data is binary.

Performance has a function called binned_residuals which groups data points by their probabilities and plots the average residual for each one.

car_glm <- glm(vs ~ wt + mpg, data = mtcars, family = "binomial") #vs is binary

plot(car_glm, which = 1) #Unhelpful

Binned residuals can help you check for randomness. I use this to visualize polynomial transformations.

library(performance)

# Fit a logistic regression model
model <- glm(as.factor(vs) ~ wt, data = mtcars, family = "binomial")

# Generate binned residuals
result <- binned_residuals(model)

# View text summary 
result
## Warning: Probably bad model fit. Only about 67% of the residuals are inside the error bounds.
# Plot the binned residuals
plot(result)