Calculating summary statistics

We have already covered this in some detail in session 3, on grouped data, so some of this code should be familiar.

We already know how to calculate summary statistics like this:

studyhabits %>% 
  summarise(mean(work_hours), mean(work_consistently), mean(grade))

  mean(work_hours) mean(work_consistently) mean(grade)
1         24.87333                3.333333    59.20667

In the workshop I showed a short-cut technique to summarise multiple variables at once. This isn’t strictly necessary, but can be a time saver.

studyhabits %>% 
  select(-unique_id) %>% 
  summarise_all(list(Mean=mean, SD=sd)) %>% 
  pivot_longer(everything()) %>% 
  arrange(name)

# A tibble: 16 × 2
   name                    value
   <chr>                   <dbl>
 1 female_Mean             0.827
 2 female_SD               0.379
 3 focus_deadline_Mean     2.85 
 4 focus_deadline_SD       0.876
 5 grade_Mean             59.2  
 6 grade_SD                9.57 
 7 msc_student_Mean        0.47 
 8 msc_student_SD          0.500
 9 progress_everyday_Mean  2.79 
10 progress_everyday_SD    0.805
11 revision_before_Mean    3.36 
12 revision_before_SD      0.853
13 work_consistently_Mean  3.33 
14 work_consistently_SD    0.851
15 work_hours_Mean        24.9  
16 work_hours_SD           4.90 

Explanation of the code We can break this down, line by line:

• First we take the studyhabits data, and remove the unique_id column using the select function.
• Next we use the summarise_all function to calculate statistics on all of these variables in one go.
• The part which says list(Mean=mean, SD=sd) tells R which functions to use. We could add the median by writing list(Mean=mean, Median=median, SD=sd)
• pivot_longer(everything()) reshapes the result from a very thin/wide table to a long-form one
• arrange(name) sorts the result

Another trick which can be useful is counting the number of responses to each option of a question. For example:

studyhabits %>% 
  group_by(work_consistently) %>% 
  count()

# A tibble: 5 × 2
# Groups:   work_consistently [5]
  work_consistently     n
              <dbl> <int>
1                 1     4
2                 2    37
3                 3   140
4                 4    93
5                 5    26

Explanation of the code: We group_by the different ‘levels’ or options of the work_consistently question. The options ranged from 1-5, so the count() function then tells us how many of each response there were.

studyhabits %>% 
  group_by(female, work_consistently) %>% 
  count()

# A tibble: 10 × 3
# Groups:   female, work_consistently [10]
   female work_consistently     n
   <lgl>              <dbl> <int>
 1 FALSE                  1     1
 2 FALSE                  2     7
 3 FALSE                  3    16
 4 FALSE                  4    19
 5 FALSE                  5     9
 6 TRUE                   1     3
 7 TRUE                   2    30
 8 TRUE                   3   124
 9 TRUE                   4    74
10 TRUE                   5    17

Confidence intervals

Once we have run a regression model, we may want confidence intervals for each coefficient. We can obtain these with tidy() in the broom package:

library(broom)

m1 <- lm(grade ~ work_hours, data=studyhabits)
m1 %>% tidy(conf.int=TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)   38.4       2.59      14.8  1.55e-37   33.3       43.5 
2 work_hours     0.838     0.102      8.19 7.68e-15    0.636      1.04

Explanation of the output:

• We add the conf.int=TRUE option to include confidence intervals
• The estimate column is the coefficient value, as you would see in the normal lm output
• The columns conf.low and conf.high represent the lower and upper bounds of the confidence interval for each coefficient.
• Ignore the other columns (p.value etc).

Standardized coefficients

As discussed in the introduction, standardizing coefficients can be helpful to enable us to compare:

• Different predictors in the same model, and
• The same predictors in different models (e.g. when the outcome measure changes)

The easiest way to standardize the coefficients from the model is to use the arm::standardize() function (in the arm package).

First we fit the model as usual:

m1 <- lm(grade ~ work_hours, data=studyhabits)

Then we have 2 options:

1. Standardize the predictor variables only
2. Standardize the predictors AND the outcome variable

m1 %>% arm::standardize()

Call:
lm(formula = grade ~ z.work_hours, data = studyhabits)

Coefficients:
 (Intercept)  z.work_hours  
      59.207         8.208  

m1 %>% arm::standardize(standardize.y=TRUE)

Call:
lm(formula = z.grade ~ z.work_hours, data = studyhabits)

Coefficients:
 (Intercept)  z.work_hours  
  -1.281e-16     4.287e-01  

Printing numbers nicely (avoiding scientific notation)

In calculating the standardized coefficients, R reported the answers using scientific notation. This can be annoying.

To fix this we can use the scipen option:

options(scipen=20)
m1 %>% 
  arm::standardize(standardize.y=TRUE) 

Call:
lm(formula = z.grade ~ z.work_hours, data = studyhabits)

Coefficients:
           (Intercept)            z.work_hours  
-0.0000000000000001281   0.4286972460463128742  

Explanation of the option name: scipen is short for scientific notation penalty. It only allows R to use scientific notation if a number is really very small, or very large. If we set scipen=20, that means a number would have to have 20 significant digits (either before or after the decimal place) before R would use the e+20 type notation.

This change persists, so you may also want to set R back to normal afterwards:

options(scipen=0)