Multivariate Statistics, Module 2: Review of Regression and ANOVA

Load in relevant packages

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.3     ✓ purrr   0.3.4
✓ tibble  3.0.5     ✓ dplyr   1.0.3
✓ tidyr   1.1.2     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.0
── Conflicts ───────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(Hmisc)

Loading required package: lattice
Loading required package: survival
Loading required package: Formula
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Attaching package: ‘Hmisc’

The following objects are masked from ‘package:dplyr’:

    src, summarize

The following objects are masked from ‘package:base’:

    format.pval, units

library(haven)

Review of ANOVA

Load in the dataset

systolic <- read_dta("systolic.dta")
View(systolic)

glimpse allows you to look more into the structure of the data:

glimpse(systolic)

Rows: 58
Columns: 3
$ drug     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ disease  <dbl> 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, …
$ systolic <dbl> 42, 44, 36, 13, 19, 22, 33, 26, 33, 21, 31, -3, 25, 25, 24, 28, 23, 34, 42, 13, 34, …

Data Cleaning and Coding:

mutate from the dplyr package in the tidyverse is a SUPER handy function that I use all the time. It allows us to generate a new variable based on a transformation (or mutation, think X-Men!) of an existing one. In this case, I am taking the original (renamed) variable drug, and I am creating a new variable called drug.fac that is equal to the old variable, but converted to a factor (that is what the as_factor function does, another one found in the tidyverse). So now, R recognizes the difference between 1 and 2 for drug.fac as actually being two different levels of a categorical variable that I can use in frequency tables, regression models, etc. I will do the same thing with the variable disease, and create a new variable called disease.fac.

Also notice that I create a whole new dataset called systolic.clean that contains these new variables. This is the version of the dataset I will use going forward. This is a really good reproducibility strategy - you don’t actually change the original raw data.

systolic.clean <- systolic %>%
  mutate(.,
         drug.fac = as_factor(drug),
         disease.fac = as_factor(disease))

describe from the Hmisc package is a nice, quick way to get descriptive statistics AND tabulations:

describe(systolic.clean)

systolic.clean 

 5  Variables      58  Observations
--------------------------------------------------------------------------------------------------------
drug : Drug Used  Format:%8.0g 
       n  missing distinct     Info     Mean      Gmd 
      58        0        4    0.936      2.5    1.305 
                                  
Value          1     2     3     4
Frequency     15    15    12    16
Proportion 0.259 0.259 0.207 0.276
--------------------------------------------------------------------------------------------------------
disease : Patient's Disease  Format:%8.0g 
       n  missing distinct     Info     Mean      Gmd 
      58        0        3    0.889    2.017    0.908 
                            
Value          1     2     3
Frequency     19    19    20
Proportion 0.328 0.328 0.345
--------------------------------------------------------------------------------------------------------
systolic : Increment in Systolic B.P.  Format:%8.0g 
       n  missing distinct     Info     Mean      Gmd      .05      .10      .25      .50      .75 
      58        0       32    0.999    18.88     14.8    -2.15     1.00     9.00    21.00    28.00 
     .90      .95 
   34.00    36.90 

lowest : -6 -5 -3 -2  1, highest: 33 34 36 42 44
--------------------------------------------------------------------------------------------------------
drug.fac 
       n  missing distinct 
      58        0        4 
                                  
Value          1     2     3     4
Frequency     15    15    12    16
Proportion 0.259 0.259 0.207 0.276
--------------------------------------------------------------------------------------------------------
disease.fac 
       n  missing distinct 
      58        0        3 
                            
Value          1     2     3
Frequency     19    19    20
Proportion 0.328 0.328 0.345
--------------------------------------------------------------------------------------------------------

Boxplots are a useful visualization tool for ANOVAs

ggplot(data = systolic.clean, mapping = aes(x = drug.fac, y = systolic)) +
  geom_boxplot() + 
  labs(title = "Distribution of Change in Systolic Pressure",
       y = "Change in Systolic Blood Pressure",
       x = "Drug Used",
       caption = "N = 58.")

Here is your basic one way ANOVA command:

anova1 <- aov(systolic ~ drug.fac, data=systolic.clean)
summary(anova1)

            Df Sum Sq Mean Sq F value   Pr(>F)    
drug.fac     3   3133  1044.4   9.086 5.75e-05 ***
Residuals   54   6207   114.9                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

This command then breaks down the mean by each drug group…

systolic.clean %>%
  group_by(drug) %>%
  summarise(n = n(), mean = mean(systolic))

Post-Hoc Tests!

Bonferroni Method

pairwise.t.test(systolic.clean$systolic, systolic.clean$drug, p.adjust.method = "bonf")

    Pairwise comparisons using t tests with pooled SD 

data:  systolic.clean$systolic and systolic.clean$drug 

  1       2       3      
2 1.00000 -       -      
3 0.00067 0.00102 -      
4 0.01154 0.01726 1.00000

P value adjustment method: bonferroni 

Tukey Method

TukeyHSD(anova1)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = systolic ~ drug.fac, data = systolic.clean)

$drug.fac
           diff        lwr       upr     p adj
2-1  -0.5333333 -10.911001  9.844334 0.9990860
3-1 -17.3166667 -28.323846 -6.309488 0.0006259
4-1 -12.5666667 -22.780896 -2.352437 0.0100845
3-2 -16.7833333 -27.790512 -5.776154 0.0009486
4-2 -12.0333333 -22.247563 -1.819104 0.0148153
4-3   4.7500000  -6.103225 15.603225 0.6542269

Calculate Effect Sizes Using the effectsize Package

Eta Squared

library(effectsize)
eta_squared(anova1, partial = FALSE)

Parameter | Eta2 |       90% CI
-------------------------------
drug.fac  | 0.34 | [0.15, 0.47]

Omega Squared

omega_squared(anova1, partial = FALSE)

Parameter | Omega2 |       90% CI
---------------------------------
drug.fac  |   0.29 | [0.11, 0.43]

Two Way ANOVA - Add in the Disease Comparison

ggplot(data = systolic.clean, mapping = aes(x = drug.fac, y = systolic, color = disease.fac)) +
  geom_boxplot() + 
  labs(title = "Distribution of Change in Systolic Pressure, by Drug Used and Disease",
       y = "Change in Systolic Blood Pressure",
       x = "Drug Used",
       caption = "N = 58.")

Main Command, with Interaction Effect

anova2 <- aov(systolic ~ drug.fac + disease.fac + drug.fac*disease.fac, data=systolic.clean)
summary(anova2)

                     Df Sum Sq Mean Sq F value   Pr(>F)    
drug.fac              3   3133  1044.4   9.456 5.58e-05 ***
disease.fac           2    419   209.4   1.896    0.162    
drug.fac:disease.fac  6    707   117.9   1.067    0.396    
Residuals            46   5081   110.5                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Effect Size

omega_squared(anova2, partial = TRUE)

Parameter            | Omega2 (partial) |       90% CI
------------------------------------------------------
drug.fac             |             0.30 | [0.11, 0.45]
disease.fac          |             0.03 | [0.00, 0.12]
drug.fac:disease.fac |         6.91e-03 | [0.00, 0.00]

Tukey HSD for Post-Hoc Tests

TukeyHSD(anova2)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = systolic ~ drug.fac + disease.fac + drug.fac * disease.fac, data = systolic.clean)

$drug.fac
           diff        lwr       upr     p adj
2-1  -0.5333333 -10.762383  9.695716 0.9990278
3-1 -17.3166667 -28.166212 -6.467121 0.0005728
4-1 -12.5666667 -22.634619 -2.498714 0.0090727
3-2 -16.7833333 -27.632879 -5.933788 0.0008634
4-2 -12.0333333 -22.101286 -1.965381 0.0133495
4-3   4.7500000  -5.947796 15.447796 0.6401821

$disease.fac
         diff       lwr      upr     p adj
2-1 -2.122807 -10.38070 6.135090 0.8085418
3-1 -6.406053 -14.56007 1.747967 0.1494034
3-2 -4.283246 -12.43727 3.870774 0.4178989

$`drug.fac:disease.fac`
               diff        lwr        upr     p adj
2:1-1:1  -1.3333333 -23.231992 20.5653251 1.0000000
3:1-1:1 -13.0000000 -38.572124 12.5721237 0.8355971
4:1-1:1 -15.7333333 -37.631992  6.1653251 0.3819307
1:2-1:1  -1.0833333 -24.427382 22.2607150 1.0000000
2:2-1:1   4.1666667 -19.177382 27.5107150 0.9999696
3:2-1:1 -24.9333333 -46.831992 -3.0346749 0.0137983
4:2-1:1 -16.5000000 -37.379552  4.3795516 0.2496012
1:3-1:1  -8.9333333 -30.831992 12.9653251 0.9568879
2:3-1:1 -11.1666667 -32.046218  9.7128849 0.7875711
3:3-1:1 -20.8333333 -44.177382  2.5107150 0.1211520
4:3-1:1 -15.1333333 -37.031992  6.7653251 0.4403630
3:1-2:1 -11.6666667 -38.077442 14.7441091 0.9271634
4:1-2:1 -14.4000000 -37.272403  8.4724028 0.5805747
1:2-2:1   0.2500000 -24.009847 24.5098466 1.0000000
2:2-2:1   5.5000000 -18.759847 29.7598466 0.9996836
3:2-2:1 -23.6000000 -46.472403 -0.7275972 0.0376995
4:2-2:1 -15.1666667 -37.065325  6.7319918 0.4370264
1:3-2:1  -7.6000000 -30.472403 15.2724028 0.9906936
2:3-2:1  -9.8333333 -31.731992 12.0653251 0.9193001
3:3-2:1 -19.5000000 -43.759847  4.7598466 0.2285024
4:3-2:1 -13.8000000 -36.672403  9.0724028 0.6411329
4:1-3:1  -2.7333333 -29.144109 23.6774424 0.9999999
1:2-3:1  11.9166667 -15.704384 39.5377171 0.9373706
2:2-3:1  17.1666667 -10.454384 44.7877171 0.5992994
3:2-3:1 -11.9333333 -38.344109 14.4774424 0.9161453
4:2-3:1  -3.5000000 -29.072124 22.0721237 0.9999980
1:3-3:1   4.0666667 -22.344109 30.4774424 0.9999932
2:3-3:1   1.8333333 -23.738790 27.4054570 1.0000000
3:3-3:1  -7.8333333 -35.454384 19.7877171 0.9975445
4:3-3:1  -2.1333333 -28.544109 24.2774424 1.0000000
1:2-4:1  14.6500000  -9.609847 38.9098466 0.6399179
2:2-4:1  19.9000000  -4.359847 44.1598466 0.2045559
3:2-4:1  -9.2000000 -32.072403 13.6724028 0.9608830
4:2-4:1  -0.7666667 -22.665325 21.1319918 1.0000000
1:3-4:1   6.8000000 -16.072403 29.6724028 0.9963056
2:3-4:1   4.5666667 -17.331992 26.4653251 0.9998583
3:3-4:1  -5.1000000 -29.359847 19.1598466 0.9998468
4:3-4:1   0.6000000 -22.272403 23.4724028 1.0000000
2:2-1:2   5.2500000 -20.322124 30.8221237 0.9998783
3:2-1:2 -23.8500000 -48.109847  0.4098466 0.0578821
4:2-1:2 -15.4166667 -38.760715  7.9273816 0.5092766
1:3-1:2  -7.8500000 -32.109847 16.4098466 0.9924943
2:3-1:2 -10.0833333 -33.427382 13.2607150 0.9368894
3:3-1:2 -19.7500000 -45.322124  5.8221237 0.2794645
4:3-1:2 -14.0500000 -38.309847 10.2098466 0.6955374
3:2-2:2 -29.1000000 -53.359847 -4.8401534 0.0075147
4:2-2:2 -20.6666667 -44.010715  2.6773816 0.1279415
1:3-2:2 -13.1000000 -37.359847 11.1598466 0.7775056
2:3-2:2 -15.3333333 -38.677382  8.0107150 0.5174698
3:3-2:2 -25.0000000 -50.572124  0.5721237 0.0606627
4:3-2:2 -19.3000000 -43.559847  4.9598466 0.2411863
4:2-3:2   8.4333333 -13.465325 30.3319918 0.9712576
1:3-3:2  16.0000000  -6.872403 38.8724028 0.4219204
2:3-3:2  13.7666667  -8.131992 35.6653251 0.5827247
3:3-3:2   4.1000000 -20.159847 28.3598466 0.9999825
4:3-3:2   9.8000000 -13.072403 32.6724028 0.9401219
1:3-4:2   7.5666667 -14.331992 29.4653251 0.9873047
2:3-4:2   5.3333333 -15.546218 26.2128849 0.9990373
3:3-4:2  -4.3333333 -27.677382 19.0107150 0.9999550
4:3-4:2   1.3666667 -20.531992 23.2653251 1.0000000
2:3-1:3  -2.2333333 -24.131992 19.6653251 0.9999999
3:3-1:3 -11.9000000 -36.159847 12.3598466 0.8644308
4:3-1:3  -6.2000000 -29.072403 16.6724028 0.9983533
3:3-2:3  -9.6666667 -33.010715 13.6773816 0.9522613
4:3-2:3  -3.9666667 -25.865325 17.9319918 0.9999647
4:3-3:3   5.7000000 -18.559847 29.9598466 0.9995565

Review of Multiple Linear Regression

New Dataset - NLSW88

nlsw88 <- read_dta("nlsw88.dta")
View(nlsw88)

glimpse(nlsw88)

Rows: 2,246
Columns: 17
$ idcode        <dbl> 1, 2, 3, 4, 6, 7, 9, 12, 13, 14, 15, 16, 18, 19, 20, 22, 23, 24, 25, 36, 39, 44…
$ age           <dbl> 37, 37, 42, 43, 42, 39, 37, 40, 40, 40, 39, 40, 40, 40, 39, 41, 42, 41, 42, 37,…
$ race          <dbl+lbl> 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ married       <dbl+lbl> 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, …
$ never_married <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0…
$ grade         <dbl> 12, 12, 12, 17, 12, 12, 12, 18, 14, 15, 16, 15, 15, 15, 15, 15, 15, 14, 14, 12,…
$ collgrad      <dbl+lbl> 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, …
$ south         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ smsa          <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, …
$ c_city        <dbl> 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1…
$ industry      <dbl+lbl>  5,  4,  4, 11,  4, 11,  5, 11, 11, 11, 11, 11,  6, 11, 11, 11, 11, 11, 11,…
$ occupation    <dbl+lbl>  6,  5,  3, 13,  6,  3,  2,  2,  3,  1,  1,  1,  5,  1,  1,  1,  1,  1,  1,…
$ union         <dbl+lbl>  1,  1, NA,  1,  0,  0,  1,  0,  0,  0,  1,  0,  0,  1,  1,  0,  0,  0,  0,…
$ wage          <dbl> 11.739125, 6.400963, 5.016723, 9.033813, 8.083731, 4.629630, 10.491142, 17.2061…
$ hours         <dbl> 48, 40, 40, 42, 48, 30, 40, 45, 8, 50, 16, 40, 40, 40, 4, 32, 45, 24, 40, 40, 4…
$ ttl_exp       <dbl> 10.333334, 13.621795, 17.730770, 13.211537, 17.820513, 7.326923, 19.044870, 15.…
$ tenure        <dbl> 5.3333335, 5.2500000, 1.2500000, 1.7500000, 17.7500000, 2.2500000, 19.0000000, …

Quick Data Cleaning - Create Factor Variable for Married

nlsw.clean <- nlsw88 %>%
  mutate(.,
         married.fac = as_factor(married))

glimpse(nlsw.clean)

Rows: 2,246
Columns: 18
$ idcode        <dbl> 1, 2, 3, 4, 6, 7, 9, 12, 13, 14, 15, 16, 18, 19, 20, 22, 23, 24, 25, 36, 39, 44…
$ age           <dbl> 37, 37, 42, 43, 42, 39, 37, 40, 40, 40, 39, 40, 40, 40, 39, 41, 42, 41, 42, 37,…
$ race          <dbl+lbl> 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ married       <dbl+lbl> 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, …
$ never_married <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0…
$ grade         <dbl> 12, 12, 12, 17, 12, 12, 12, 18, 14, 15, 16, 15, 15, 15, 15, 15, 15, 14, 14, 12,…
$ collgrad      <dbl+lbl> 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, …
$ south         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ smsa          <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, …
$ c_city        <dbl> 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1…
$ industry      <dbl+lbl>  5,  4,  4, 11,  4, 11,  5, 11, 11, 11, 11, 11,  6, 11, 11, 11, 11, 11, 11,…
$ occupation    <dbl+lbl>  6,  5,  3, 13,  6,  3,  2,  2,  3,  1,  1,  1,  5,  1,  1,  1,  1,  1,  1,…
$ union         <dbl+lbl>  1,  1, NA,  1,  0,  0,  1,  0,  0,  0,  1,  0,  0,  1,  1,  0,  0,  0,  0,…
$ wage          <dbl> 11.739125, 6.400963, 5.016723, 9.033813, 8.083731, 4.629630, 10.491142, 17.2061…
$ hours         <dbl> 48, 40, 40, 42, 48, 30, 40, 45, 8, 50, 16, 40, 40, 40, 4, 32, 45, 24, 40, 40, 4…
$ ttl_exp       <dbl> 10.333334, 13.621795, 17.730770, 13.211537, 17.820513, 7.326923, 19.044870, 15.…
$ tenure        <dbl> 5.3333335, 5.2500000, 1.2500000, 1.7500000, 17.7500000, 2.2500000, 19.0000000, …
$ married.fac   <fct> single, single, single, married, married, married, single, married, married, ma…

Scatterplot of Wage and Tenure

ggplot(data = nlsw.clean, mapping = aes(x = tenure, y = wage)) + 
  geom_point(colour = "blue", size = 1.5) +
  labs(title = "Relationship between Wage and Tenure",
       x = "Years in Job", 
       y = "Hourly Wage",
       caption = "Data from the National Longitudinal Survey of Women (1988). N = 2,246.")

Simple (Bivariate) Linear Regression

linearmodel1 <- lm(wage ~ tenure, data = nlsw.clean)
summary(linearmodel1)

Call:
lm(formula = wage ~ tenure, data = nlsw.clean)

Residuals:
   Min     1Q Median     3Q    Max 
-8.027 -3.187 -1.470  1.448 33.786 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.68132    0.17726  37.692   <2e-16 ***
tenure       0.18587    0.02181   8.524   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.674 on 2229 degrees of freedom
  (15 observations deleted due to missingness)
Multiple R-squared:  0.03157,   Adjusted R-squared:  0.03114 
F-statistic: 72.66 on 1 and 2229 DF,  p-value: < 2.2e-16

Scatterplot With Added Regression Line (stat_smooth)

ggplot(data = nlsw.clean, mapping = aes(x = tenure, y = wage)) + 
  geom_point(colour = "blue", size = 1.5) +
  labs(title = "Relationship between Wage and Tenure",
       x = "Years in Job", 
       y = "Hourly Wage",
       caption = "Data from the National Longitudinal Survey of Women (1988). N = 2,246. Red line represents a smoothed regression line.") + 
  stat_smooth(method="loess", col="red", size=1)

Adding More Predictors- Now We Have a Multiple Linear Regression (MLR)

linearmodel2 <- lm(wage ~ tenure + age + grade + married.fac, data = nlsw.clean)
summary(linearmodel2)

Call:
lm(formula = wage ~ tenure + age + grade + married.fac, data = nlsw.clean)

Residuals:
    Min      1Q  Median      3Q     Max 
-10.921  -2.726  -1.084   1.040  35.362 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         0.78784    1.60824   0.490   0.6243    
tenure              0.14849    0.02097   7.082  1.9e-12 ***
age                -0.07066    0.03742  -1.888   0.0591 .  
grade               0.70217    0.04561  15.393  < 2e-16 ***
married.facmarried -0.48747    0.23828  -2.046   0.0409 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.388 on 2224 degrees of freedom
  (17 observations deleted due to missingness)
Multiple R-squared:  0.1284,    Adjusted R-squared:  0.1268 
F-statistic: 81.92 on 4 and 2224 DF,  p-value: < 2.2e-16

Need Standardized (Beta) Coefficients? Try the reghelper Package…

library(reghelper)
lm2.beta <- beta(linearmodel2)
lm2.beta

Call:
lm(formula = "wage.z ~ tenure.z + age.z + grade.z + married.facmarried.z", 
    data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.8938 -0.4726 -0.1879  0.1804  6.1323 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)          -2.157e-17  1.979e-02   0.000   1.0000    
tenure.z              1.418e-01  2.002e-02   7.082  1.9e-12 ***
age.z                -3.752e-02  1.987e-02  -1.888   0.0591 .  
grade.z               3.074e-01  1.997e-02  15.393  < 2e-16 ***
married.facmarried.z -4.054e-02  1.981e-02  -2.046   0.0409 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9344 on 2224 degrees of freedom
Multiple R-squared:  0.1284,    Adjusted R-squared:  0.1268 
F-statistic: 81.92 on 4 and 2224 DF,  p-value: < 2.2e-16

Need a Nice Table to Summarize? Try the modelsummary Package

library(modelsummary)
models <- list(linearmodel1, linearmodel2)
modelsummary(models, output = "default")

	Model 1	Model 2
(Intercept)	6.681	0.788
	(0.177)	(1.608)
tenure	0.186	0.148
	(0.022)	(0.021)
age		-0.071
		(0.037)
grade		0.702
		(0.046)
married.facmarried		-0.487
		(0.238)
Num.Obs.	2231	2229
R2	0.032	0.128
R2 Adj.	0.031	0.127
AIC	14080.9	13841.0
BIC	14098.0	13875.3
Log.Lik.	-7037.454	-6914.507
F	72.663	81.919

You Can Also Make a Model Plot!

modelplot(lm2.beta)