Init

library(kirkegaard)

## Loading required package: tidyverse

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6      ✔ purrr   0.3.5 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## Loading required package: magrittr
## 
## 
## Attaching package: 'magrittr'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     set_names
## 
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## 
## Loading required package: weights
## 
## Loading required package: Hmisc
## 
## Loading required package: lattice
## 
## Loading required package: survival
## 
## Loading required package: Formula
## 
## 
## Attaching package: 'Hmisc'
## 
## 
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## 
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## 
## 
## Loading required package: assertthat
## 
## 
## Attaching package: 'assertthat'
## 
## 
## The following object is masked from 'package:tibble':
## 
##     has_name
## 
## 
## Loading required package: psych
## 
## 
## Attaching package: 'psych'
## 
## 
## The following object is masked from 'package:Hmisc':
## 
##     describe
## 
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
## 
## 
## 
## Attaching package: 'kirkegaard'
## 
## 
## The following object is masked from 'package:psych':
## 
##     rescale
## 
## 
## The following object is masked from 'package:assertthat':
## 
##     are_equal
## 
## 
## The following objects are masked from 'package:purrr':
## 
##     is_logical, is_numeric
## 
## 
## The following object is masked from 'package:base':
## 
##     +

load_packages(
  rms,
  ggeffects,
  mirt
)

## Loading required package: SparseM
## 
## Attaching package: 'SparseM'
## 
## The following object is masked from 'package:base':
## 
##     backsolve
## 
## Loading required package: stats4

theme_set(theme_bw())

options(
    digits = 3
)

Functions

#compare total effect of groups of variables using ANOVA eta
compare_total_effects = function(model, groups) {
  #fit anova
  anova_fit = car::Anova(model)
  
  #omega squared
  anova_stats = sjstats::anova_stats(anova_fit)
  
  #sum effects by group
  list(
    anova_stats = anova_stats,
    sums = map2_df(groups, names(groups), function(group, group_name) {
      group2 = group
      tibble(
        group = group_name,
        sum_omegasq = anova_stats %>% filter(term %in% group2) %>% pull(omegasq) %>% sum(),
        sum_omega = sqrt(sum_omegasq)
      )
    })
  )
}

mpg$year %<>% factor()
compare_total_effects(
  lm(displ ~ year + cyl + drv + cty + hwy, data = mpg),
  list(group1 = c("year", "cyl"), group2 = c("drv", "cty", "hwy"))
)

## $anova_stats
## term      |  sumsq | meansq |  df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ------------------------------------------------------------------------------------------------------------------------------------------
## year      |  0.823 |  0.823 |   1 |     4.731 |   0.031 | 0.007 |         0.020 |   0.006 |           0.016 |     0.006 |    0.144 | 0.585
## cyl       | 61.009 | 61.009 |   1 |   350.607 |  < .001 | 0.551 |         0.607 |   0.549 |           0.599 |     0.550 |    1.243 | 1.000
## drv       |  8.928 |  4.464 |   2 |    25.654 |  < .001 | 0.081 |         0.184 |   0.077 |           0.174 |     0.078 |    0.475 | 1.000
## cty       |  0.003 |  0.003 |   1 |     0.016 |   0.898 | 0.000 |         0.000 |  -0.002 |          -0.004 |    -0.002 |    0.009 | 0.052
## hwy       |  0.402 |  0.402 |   1 |     2.308 |   0.130 | 0.004 |         0.010 |   0.002 |           0.006 |     0.002 |    0.101 | 0.330
## Residuals | 39.500 |  0.174 | 227 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 2 × 3
##   group  sum_omegasq sum_omega
##   <chr>        <dbl>     <dbl>
## 1 group1       0.555     0.745
## 2 group2       0.077     0.277

Data

d = read_rds("data/data_out.rds")
d_table = df_var_table(d)

aldo = readxl::read_excel("data/Aldo_dysgenics_VES.xlsx")
aldo_table = df_var_table(aldo)

d$SOCSTAT = aldo$SOCSTAT
d$SOCSTAT2 = aldo$SOCSTAT2

Recode

#intelligence g
#g tests
g_tests_early = c("VE_time1", "AR_time1", "PA", "GIT", "AFQT")
g_tests_later = c("VE_time2", "AR_time2", "WAIS_BD", "WAIS_GI", "WRAT", "PASAT", "WLGT", "copy_direct", "copy_immediate", "copy_delayed", "CVLT", "WCST", "GPT_left", "GPT_right")
g_tests = c(g_tests_early, g_tests_later)

#test correlations
wtd.cors(d[g_tests])

##                VE_time1 AR_time1    PA   GIT  AFQT VE_time2 AR_time2 WAIS_BD
## VE_time1          1.000    0.699 0.516 0.659 0.714    0.824    0.642   0.437
## AR_time1          0.699    1.000 0.576 0.589 0.737    0.658    0.785   0.502
## PA                0.516    0.576 1.000 0.467 0.728    0.484    0.545   0.634
## GIT               0.659    0.589 0.467 1.000 0.645    0.620    0.548   0.418
## AFQT              0.714    0.737 0.728 0.645 1.000    0.670    0.688   0.629
## VE_time2          0.824    0.658 0.484 0.620 0.670    1.000    0.691   0.453
## AR_time2          0.642    0.785 0.545 0.548 0.688    0.691    1.000   0.532
## WAIS_BD           0.437    0.502 0.634 0.418 0.629    0.453    0.532   1.000
## WAIS_GI           0.725    0.635 0.482 0.582 0.626    0.719    0.622   0.453
## WRAT              0.746    0.589 0.412 0.517 0.578    0.766    0.585   0.382
## PASAT             0.408    0.521 0.371 0.365 0.432    0.440    0.562   0.388
## WLGT              0.443    0.370 0.289 0.310 0.360    0.463    0.365   0.281
## copy_direct       0.290    0.333 0.380 0.254 0.372    0.325    0.384   0.398
## copy_immediate    0.303    0.335 0.456 0.311 0.452    0.316    0.387   0.490
## copy_delayed      0.301    0.337 0.452 0.308 0.451    0.315    0.385   0.492
## CVLT              0.317    0.331 0.264 0.250 0.312    0.333    0.356   0.269
## WCST              0.327    0.360 0.331 0.282 0.368    0.361    0.396   0.356
## GPT_left          0.208    0.212 0.263 0.192 0.269    0.226    0.236   0.307
## GPT_right         0.204    0.201 0.261 0.174 0.253    0.220    0.240   0.301
##                WAIS_GI  WRAT PASAT  WLGT copy_direct copy_immediate
## VE_time1         0.725 0.746 0.408 0.443       0.290          0.303
## AR_time1         0.635 0.589 0.521 0.370       0.333          0.335
## PA               0.482 0.412 0.371 0.289       0.380          0.456
## GIT              0.582 0.517 0.365 0.310       0.254          0.311
## AFQT             0.626 0.578 0.432 0.360       0.372          0.452
## VE_time2         0.719 0.766 0.440 0.463       0.325          0.316
## AR_time2         0.622 0.585 0.562 0.365       0.384          0.387
## WAIS_BD          0.453 0.382 0.388 0.281       0.398          0.490
## WAIS_GI          1.000 0.652 0.366 0.414       0.278          0.343
## WRAT             0.652 1.000 0.417 0.504       0.269          0.269
## PASAT            0.366 0.417 1.000 0.357       0.247          0.282
## WLGT             0.414 0.504 0.357 1.000       0.176          0.213
## copy_direct      0.278 0.269 0.247 0.176       1.000          0.474
## copy_immediate   0.343 0.269 0.282 0.213       0.474          1.000
## copy_delayed     0.345 0.271 0.281 0.216       0.482          0.915
## CVLT             0.329 0.309 0.289 0.278       0.205          0.318
## WCST             0.330 0.292 0.285 0.209       0.291          0.268
## GPT_left         0.186 0.204 0.216 0.157       0.223          0.227
## GPT_right        0.173 0.196 0.226 0.168       0.223          0.194
##                copy_delayed  CVLT  WCST GPT_left GPT_right
## VE_time1              0.301 0.317 0.327    0.208     0.204
## AR_time1              0.337 0.331 0.360    0.212     0.201
## PA                    0.452 0.264 0.331    0.263     0.261
## GIT                   0.308 0.250 0.282    0.192     0.174
## AFQT                  0.451 0.312 0.368    0.269     0.253
## VE_time2              0.315 0.333 0.361    0.226     0.220
## AR_time2              0.385 0.356 0.396    0.236     0.240
## WAIS_BD               0.492 0.269 0.356    0.307     0.301
## WAIS_GI               0.345 0.329 0.330    0.186     0.173
## WRAT                  0.271 0.309 0.292    0.204     0.196
## PASAT                 0.281 0.289 0.285    0.216     0.226
## WLGT                  0.216 0.278 0.209    0.157     0.168
## copy_direct           0.482 0.205 0.291    0.223     0.223
## copy_immediate        0.915 0.318 0.268    0.227     0.194
## copy_delayed          1.000 0.327 0.268    0.227     0.199
## CVLT                  0.327 1.000 0.192    0.115     0.117
## WCST                  0.268 0.192 1.000    0.196     0.186
## GPT_left              0.227 0.115 0.196    1.000     0.634
## GPT_right             0.199 0.117 0.186    0.634     1.000

#impute missing data
g_all_tests = d[g_tests] %>% miss_impute()

#all tests g
fa_g = fa(g_all_tests)
fa_g

## Factor Analysis using method =  minres
## Call: fa(r = g_all_tests)
## Standardized loadings (pattern matrix) based upon correlation matrix
##                 MR1   h2   u2 com
## VE_time1       0.82 0.67 0.33   1
## AR_time1       0.81 0.66 0.34   1
## PA             0.71 0.50 0.50   1
## GIT            0.69 0.48 0.52   1
## AFQT           0.85 0.73 0.27   1
## VE_time2       0.82 0.67 0.33   1
## AR_time2       0.82 0.67 0.33   1
## WAIS_BD        0.67 0.45 0.55   1
## WAIS_GI        0.76 0.58 0.42   1
## WRAT           0.73 0.53 0.47   1
## PASAT          0.57 0.32 0.68   1
## WLGT           0.49 0.24 0.76   1
## copy_direct    0.47 0.22 0.78   1
## copy_immediate 0.55 0.30 0.70   1
## copy_delayed   0.55 0.30 0.70   1
## CVLT           0.42 0.18 0.82   1
## WCST           0.46 0.21 0.79   1
## GPT_left       0.34 0.12 0.88   1
## GPT_right      0.33 0.11 0.89   1
## 
##                 MR1
## SS loadings    7.94
## Proportion Var 0.42
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 factor is sufficient.
## 
## The degrees of freedom for the null model are  171  and the objective function was  12.6 with Chi Square of  56117
## The degrees of freedom for the model are 152  and the objective function was  3.89 
## 
## The root mean square of the residuals (RMSR) is  0.09 
## The df corrected root mean square of the residuals is  0.1 
## 
## The harmonic number of observations is  4462 with the empirical chi square  13120  with prob <  0 
## The total number of observations was  4462  with Likelihood Chi Square =  17324  with prob <  0 
## 
## Tucker Lewis Index of factoring reliability =  0.655
## RMSEA index =  0.159  and the 90 % confidence intervals are  0.157 0.161
## BIC =  16047
## Fit based upon off diagonal values = 0.95
## Measures of factor score adequacy             
##                                                    MR1
## Correlation of (regression) scores with factors   0.97
## Multiple R square of scores with factors          0.95
## Minimum correlation of possible factor scores     0.89

fa_g$loadings[, 1] %>% describe() %>% as.matrix()

##    vars  n  mean    sd median trimmed   mad  min   max range   skew kurtosis
## X1    1 19 0.624 0.172  0.673   0.628 0.213 0.33 0.853 0.524 -0.238    -1.44
##        se
## X1 0.0395

#earlier tests only
#impute missing data
g_early_tests = d[g_tests_early] %>% miss_impute()
fa_g_time1 = fa(g_early_tests)
fa_g_time1

## Factor Analysis using method =  minres
## Call: fa(r = g_early_tests)
## Standardized loadings (pattern matrix) based upon correlation matrix
##           MR1   h2   u2 com
## VE_time1 0.82 0.67 0.33   1
## AR_time1 0.82 0.68 0.32   1
## PA       0.70 0.49 0.51   1
## GIT      0.73 0.53 0.47   1
## AFQT     0.92 0.85 0.15   1
## 
##                 MR1
## SS loadings    3.21
## Proportion Var 0.64
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 factor is sufficient.
## 
## The degrees of freedom for the null model are  10  and the objective function was  3.11 with Chi Square of  13878
## The degrees of freedom for the model are 5  and the objective function was  0.16 
## 
## The root mean square of the residuals (RMSR) is  0.04 
## The df corrected root mean square of the residuals is  0.06 
## 
## The harmonic number of observations is  4388 with the empirical chi square  170  with prob <  6.6e-35 
## The total number of observations was  4462  with Likelihood Chi Square =  697  with prob <  2.6e-148 
## 
## Tucker Lewis Index of factoring reliability =  0.9
## RMSEA index =  0.176  and the 90 % confidence intervals are  0.165 0.187
## BIC =  655
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    MR1
## Correlation of (regression) scores with factors   0.96
## Multiple R square of scores with factors          0.93
## Minimum correlation of possible factor scores     0.85

fa_g_time1$loadings[, 1] %>% describe() %>% as.matrix()

##    vars n  mean     sd median trimmed   mad   min  max range  skew kurtosis
## X1    1 5 0.797 0.0873  0.816   0.797 0.134 0.701 0.92 0.219 0.193    -1.83
##       se
## X1 0.039

#save scores, standardize to white
d$g = fa_g$scores[, 1] %>% standardize(focal_group = d$race == "White")
d$g_time1 = fa_g_time1$scores[, 1] %>% standardize(focal_group = d$race == "White")

#standard MMPI scales
d %<>% mutate(
  MMPI_L = MM010568,
  MMPI_F = MM010569,
  MMPI_K = MM010570,
  MMPI_HS = MM010571,
  MMPI_D = MM010572,
  MMPI_HY = MM010573,
  MMPI_PD = MM010574,
  MMPI_MF = MM010575,
  MMPI_PA = MM010576,
  MMPI_PT = MM010577,
  MMPI_SC = MM010578,
  MMPI_MA = MM010579,
  MMPI_SI = MM010580,
  MMPI_ES = MM010581,
  
  #g
  g = standardize(g, focal_group = d$race=="White"),
  
  #standardize outcomes
  income = income %>% standardize(focal_group = d$race=="White"),
  education = education %>% standardize(focal_group = d$race=="White"),
  occu_status = OCCSTAT %>% standardize(focal_group = d$race=="White"),
  unemployment = unemployment_3yrs %>% standardize(focal_group = d$race=="White"),
  discharge_ok = (DISCHG2 %in% c(1, 2)),
  married = (MARISTAT == 1),
  military_rank = RANK %>% standardize(focal_group = d$race=="White")
)

#MMPI vars
MMPI_vars = d %>% select(MMPI_L:MMPI_ES) %>% names()

#standardize to white norms
for (v in MMPI_vars) d[[v]] %<>% standardize(focal_group = d$race=="White")

#MMPI general factor
mmpi_irt = mirt(
  d %>% select(MM010001:MM010566),
  model = 1,
  TOL = .0005,
  verbose = F
)

#scores
mmpi_irt_scores = fscores(mmpi_irt, full.scores = T, full.scores.SE = T)
d$MMPI_p = standardize(mmpi_irt_scores[, 1], focal_group = d$race == "White")

Analysis

MMPI scales vs. g

#predict income
list(
  ols(income ~ g, data = d),
  ols(as.formula(str_glue("income ~ {str_c(MMPI_vars, collapse = ' + ')}")), data = d),
  ols(as.formula(str_glue("income ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), data = d)
  ) %>% 
  summarize_models()

#education
list(
  ols(education ~ g, data = d),
  ols(as.formula(str_glue("education ~ {str_c(MMPI_vars, collapse = ' + ')}")), data = d),
  ols(as.formula(str_glue("education ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), data = d)
  ) %>% 
  summarize_models()

#occupational status
# list(
#   ols(occu_status ~ g, data = d),
#   ols(as.formula(str_glue("occu_status ~ {str_c(MMPI_vars, collapse = ' + ')}")), data = d),
#   ols(as.formula(str_glue("occu_status ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), data = d)
#   ) %>% 
#   summarize_models()

#unemployment
list(
  ols(unemployment ~ g, data = d),
  ols(as.formula(str_glue("unemployment ~ {str_c(MMPI_vars, collapse = ' + ')}")), data = d),
  ols(as.formula(str_glue("unemployment ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), data = d)
  ) %>% 
  summarize_models()

#military rank
list(
  ols(military_rank ~ g, data = d),
  ols(as.formula(str_glue("military_rank ~ {str_c(MMPI_vars, collapse = ' + ')}")), data = d),
  ols(as.formula(str_glue("military_rank ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), data = d)
  ) %>% 
  summarize_models()

#discharge status
list(
  lrm(discharge_ok ~ g, data = d),
  lrm(as.formula(str_glue("discharge_ok ~ {str_c(MMPI_vars, collapse = ' + ')}")), data = d),
  lrm(as.formula(str_glue("discharge_ok ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), data = d)
  )

## [[1]]
## Logistic Regression Model
##  
##  lrm(formula = discharge_ok ~ g, data = d)
##  
##                         Model Likelihood      Discrimination    Rank Discrim.    
##                               Ratio Test             Indexes          Indexes    
##  Obs          4462    LR chi2      39.50      R2       0.052    C       0.701    
##   FALSE         84    d.f.             1     R2(1,4462)0.009    Dxy     0.402    
##   TRUE        4378    Pr(> chi2) <0.0001    R2(1,247.3)0.144    gamma   0.404    
##  max |deriv| 2e-14                            Brier    0.018    tau-a   0.015    
##  
##            Coef   S.E.   Wald Z Pr(>|Z|)
##  Intercept 4.3110 0.1461 29.51  <0.0001 
##  g         0.6463 0.1051  6.15  <0.0001 
##  
## 
## [[2]]
## Logistic Regression Model
##  
##  lrm(formula = as.formula(str_glue("discharge_ok ~ {str_c(MMPI_vars, collapse = ' + ')}")), 
##      data = d)
##  
##                         Model Likelihood       Discrimination    Rank Discrim.    
##                               Ratio Test              Indexes          Indexes    
##  Obs          4462    LR chi2      54.38       R2       0.071    C       0.725    
##   FALSE         84    d.f.            14     R2(14,4462)0.009    Dxy     0.450    
##   TRUE        4378    Pr(> chi2) <0.0001    R2(14,247.3)0.151    gamma   0.451    
##  max |deriv| 4e-09                             Brier    0.018    tau-a   0.017    
##  
##            Coef    S.E.   Wald Z Pr(>|Z|)
##  Intercept  4.3651 0.1501 29.08  <0.0001 
##  MMPI_L    -0.2666 0.1183 -2.25  0.0242  
##  MMPI_F    -0.1077 0.1749 -0.62  0.5380  
##  MMPI_K     0.4514 0.2142  2.11  0.0351  
##  MMPI_HS   -0.1093 0.1900 -0.58  0.5650  
##  MMPI_D    -0.0742 0.1933 -0.38  0.7012  
##  MMPI_HY    0.2758 0.1929  1.43  0.1529  
##  MMPI_PD   -0.5328 0.1501 -3.55  0.0004  
##  MMPI_MF    0.2218 0.1246  1.78  0.0750  
##  MMPI_PA    0.2514 0.1474  1.71  0.0881  
##  MMPI_PT    0.0411 0.2288  0.18  0.8576  
##  MMPI_SC    0.0319 0.2457  0.13  0.8966  
##  MMPI_MA   -0.2854 0.1540 -1.85  0.0639  
##  MMPI_SI    0.1497 0.1986  0.75  0.4510  
##  MMPI_ES    0.0391 0.1859  0.21  0.8335  
##  
## 
## [[3]]
## Logistic Regression Model
##  
##  lrm(formula = as.formula(str_glue("discharge_ok ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), 
##      data = d)
##  
##                         Model Likelihood       Discrimination    Rank Discrim.    
##                               Ratio Test              Indexes          Indexes    
##  Obs          4462    LR chi2      71.15       R2       0.093    C       0.756    
##   FALSE         84    d.f.            15     R2(15,4462)0.013    Dxy     0.511    
##   TRUE        4378    Pr(> chi2) <0.0001    R2(15,247.3)0.203    gamma   0.513    
##  max |deriv| 3e-08                             Brier    0.018    tau-a   0.019    
##  
##            Coef    S.E.   Wald Z Pr(>|Z|)
##  Intercept  4.5296 0.1655 27.36  <0.0001 
##  g          0.5596 0.1391  4.02  <0.0001 
##  MMPI_L    -0.0985 0.1256 -0.78  0.4330  
##  MMPI_F    -0.1043 0.1753 -0.60  0.5517  
##  MMPI_K     0.3838 0.2172  1.77  0.0772  
##  MMPI_HS   -0.0687 0.1911 -0.36  0.7190  
##  MMPI_D     0.0001 0.1939  0.00  0.9995  
##  MMPI_HY    0.1524 0.1961  0.78  0.4372  
##  MMPI_PD   -0.5377 0.1521 -3.53  0.0004  
##  MMPI_MF    0.0301 0.1334  0.23  0.8216  
##  MMPI_PA    0.1999 0.1479  1.35  0.1767  
##  MMPI_PT    0.0502 0.2284  0.22  0.8260  
##  MMPI_SC    0.0693 0.2474  0.28  0.7793  
##  MMPI_MA   -0.2432 0.1531 -1.59  0.1122  
##  MMPI_SI    0.0563 0.2013  0.28  0.7798  
##  MMPI_ES   -0.2050 0.1951 -1.05  0.2935  
##

#marrital status
list(
  lrm(married ~ g, data = d),
  lrm(as.formula(str_glue("married ~ {str_c(MMPI_vars, collapse = ' + ')}")), data = d),
  lrm(as.formula(str_glue("married ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), data = d)
  )

## [[1]]
## Frequencies of Missing Values Due to Each Variable
## married       g 
##       3       0 
## 
## Logistic Regression Model
##  
##  lrm(formula = married ~ g, data = d)
##  
##  
##                         Model Likelihood       Discrimination    Rank Discrim.    
##                               Ratio Test              Indexes          Indexes    
##  Obs          4459    LR chi2      44.35       R2       0.014    C       0.561    
##   FALSE       1184    d.f.             1      R2(1,4459)0.010    Dxy     0.123    
##   TRUE        3275    Pr(> chi2) <0.0001    R2(1,2608.8)0.016    gamma   0.123    
##  max |deriv| 4e-14                             Brier    0.193    tau-a   0.048    
##  
##            Coef   S.E.   Wald Z Pr(>|Z|)
##  Intercept 1.0698 0.0354 30.20  <0.0001 
##  g         0.2107 0.0318  6.64  <0.0001 
##  
## 
## [[2]]
## Frequencies of Missing Values Due to Each Variable
## married  MMPI_L  MMPI_F  MMPI_K MMPI_HS  MMPI_D MMPI_HY MMPI_PD MMPI_MF MMPI_PA 
##       3       0       0       0       0       0       0       0       0       0 
## MMPI_PT MMPI_SC MMPI_MA MMPI_SI MMPI_ES 
##       0       0       0       0       0 
## 
## Logistic Regression Model
##  
##  lrm(formula = as.formula(str_glue("married ~ {str_c(MMPI_vars, collapse = ' + ')}")), 
##      data = d)
##  
##  
##                         Model Likelihood        Discrimination    Rank Discrim.    
##                               Ratio Test               Indexes          Indexes    
##  Obs          4459    LR chi2     395.26        R2       0.124    C       0.689    
##   FALSE       1184    d.f.            14      R2(14,4459)0.082    Dxy     0.378    
##   TRUE        3275    Pr(> chi2) <0.0001    R2(14,2608.8)0.136    gamma   0.378    
##  max |deriv| 1e-08                              Brier    0.178    tau-a   0.147    
##  
##            Coef    S.E.   Wald Z Pr(>|Z|)
##  Intercept  1.1866 0.0384 30.89  <0.0001 
##  MMPI_L    -0.2369 0.0392 -6.04  <0.0001 
##  MMPI_F    -0.0276 0.0614 -0.45  0.6526  
##  MMPI_K     0.1014 0.0693  1.46  0.1437  
##  MMPI_HS    0.0583 0.0629  0.93  0.3544  
##  MMPI_D    -0.4114 0.0641 -6.42  <0.0001 
##  MMPI_HY    0.3516 0.0636  5.53  <0.0001 
##  MMPI_PD   -0.3288 0.0504 -6.53  <0.0001 
##  MMPI_MF   -0.2990 0.0394 -7.58  <0.0001 
##  MMPI_PA    0.0148 0.0488  0.30  0.7622  
##  MMPI_PT    0.2143 0.0760  2.82  0.0048  
##  MMPI_SC   -0.2545 0.0826 -3.08  0.0021  
##  MMPI_MA   -0.1484 0.0506 -2.93  0.0034  
##  MMPI_SI    0.3610 0.0645  5.60  <0.0001 
##  MMPI_ES    0.0624 0.0611  1.02  0.3070  
##  
## 
## [[3]]
## Frequencies of Missing Values Due to Each Variable
## married       g  MMPI_L  MMPI_F  MMPI_K MMPI_HS  MMPI_D MMPI_HY MMPI_PD MMPI_MF 
##       3       0       0       0       0       0       0       0       0       0 
## MMPI_PA MMPI_PT MMPI_SC MMPI_MA MMPI_SI MMPI_ES 
##       0       0       0       0       0       0 
## 
## Logistic Regression Model
##  
##  lrm(formula = as.formula(str_glue("married ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), 
##      data = d)
##  
##  
##                         Model Likelihood        Discrimination    Rank Discrim.    
##                               Ratio Test               Indexes          Indexes    
##  Obs          4459    LR chi2     421.30        R2       0.131    C       0.696    
##   FALSE       1184    d.f.            15      R2(15,4459)0.087    Dxy     0.392    
##   TRUE        3275    Pr(> chi2) <0.0001    R2(15,2608.8)0.144    gamma   0.392    
##  max |deriv| 2e-08                              Brier    0.177    tau-a   0.153    
##  
##            Coef    S.E.   Wald Z Pr(>|Z|)
##  Intercept  1.2162 0.0392 31.02  <0.0001 
##  g          0.2171 0.0427  5.08  <0.0001 
##  MMPI_L    -0.1730 0.0412 -4.20  <0.0001 
##  MMPI_F    -0.0314 0.0616 -0.51  0.6108  
##  MMPI_K     0.0637 0.0699  0.91  0.3622  
##  MMPI_HS    0.0836 0.0633  1.32  0.1869  
##  MMPI_D    -0.3929 0.0643 -6.11  <0.0001 
##  MMPI_HY    0.3077 0.0643  4.78  <0.0001 
##  MMPI_PD   -0.3242 0.0505 -6.41  <0.0001 
##  MMPI_MF   -0.3707 0.0421 -8.81  <0.0001 
##  MMPI_PA   -0.0131 0.0492 -0.27  0.7908  
##  MMPI_PT    0.2209 0.0761  2.90  0.0037  
##  MMPI_SC   -0.2400 0.0828 -2.90  0.0037  
##  MMPI_MA   -0.1336 0.0507 -2.63  0.0084  
##  MMPI_SI    0.3282 0.0650  5.05  <0.0001 
##  MMPI_ES   -0.0297 0.0639 -0.47  0.6415  
##

#omega comparisons
groups1 = list(g = "g", MMPI = MMPI_vars)
compare_total_effects(
  lm(as.formula(str_glue("income ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), data = d),
  groups1
)

## $anova_stats
## term      |    sumsq |  meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ----------------------------------------------------------------------------------------------------------------------------------------------
## g         |  255.817 | 255.817 |    1 |   305.486 |  < .001 | 0.064 |         0.065 |   0.063 |           0.065 |     0.063 |    0.265 | 1.000
## MMPI_L    |   16.898 |  16.898 |    1 |    20.179 |  < .001 | 0.004 |         0.005 |   0.004 |           0.004 |     0.004 |    0.068 | 0.994
## MMPI_F    |    2.796 |   2.796 |    1 |     3.339 |   0.068 | 0.001 |         0.001 |   0.000 |           0.001 |     0.000 |    0.028 | 0.447
## MMPI_K    |    4.566 |   4.566 |    1 |     5.452 |   0.020 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.035 | 0.646
## MMPI_HS   |    1.330 |   1.330 |    1 |     1.588 |   0.208 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.019 | 0.243
## MMPI_D    |    0.169 |   0.169 |    1 |     0.201 |   0.654 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.007 | 0.073
## MMPI_HY   |    5.970 |   5.970 |    1 |     7.129 |   0.008 | 0.001 |         0.002 |   0.001 |           0.001 |     0.001 |    0.040 | 0.761
## MMPI_PD   |   56.599 |  56.599 |    1 |    67.589 |  < .001 | 0.014 |         0.015 |   0.014 |           0.015 |     0.014 |    0.125 | 1.000
## MMPI_MF   |    1.091 |   1.091 |    1 |     1.303 |   0.254 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.017 | 0.208
## MMPI_PA   |    2.143 |   2.143 |    1 |     2.559 |   0.110 | 0.001 |         0.001 |   0.000 |           0.000 |     0.000 |    0.024 | 0.359
## MMPI_PT   |    3.901 |   3.901 |    1 |     4.658 |   0.031 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.033 | 0.579
## MMPI_SC   |    1.586 |   1.586 |    1 |     1.894 |   0.169 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.021 | 0.280
## MMPI_MA   |    3.857 |   3.857 |    1 |     4.606 |   0.032 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.033 | 0.574
## MMPI_SI   |    0.034 |   0.034 |    1 |     0.040 |   0.841 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.003 | 0.055
## MMPI_ES   |   16.158 |  16.158 |    1 |    19.296 |  < .001 | 0.004 |         0.004 |   0.004 |           0.004 |     0.004 |    0.067 | 0.993
## Residuals | 3651.117 |   0.837 | 4360 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 2 × 3
##   group sum_omegasq sum_omega
##   <chr>       <dbl>     <dbl>
## 1 g           0.063     0.251
## 2 MMPI        0.026     0.161

compare_total_effects(
  lm(as.formula(str_glue("education ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), data = d),
  groups1
)

## $anova_stats
## term      |    sumsq |  meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ----------------------------------------------------------------------------------------------------------------------------------------------
## g         |  567.508 | 567.508 |    1 |   902.796 |  < .001 | 0.161 |         0.169 |   0.160 |           0.168 |     0.161 |    0.451 | 1.000
## MMPI_L    |    0.513 |   0.513 |    1 |     0.816 |   0.366 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.014 | 0.147
## MMPI_F    |    3.410 |   3.410 |    1 |     5.425 |   0.020 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.035 | 0.644
## MMPI_K    |   20.984 |  20.984 |    1 |    33.382 |  < .001 | 0.006 |         0.007 |   0.006 |           0.007 |     0.006 |    0.087 | 1.000
## MMPI_HS   |    4.108 |   4.108 |    1 |     6.535 |   0.011 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.038 | 0.725
## MMPI_D    |    5.962 |   5.962 |    1 |     9.484 |   0.002 | 0.002 |         0.002 |   0.002 |           0.002 |     0.002 |    0.046 | 0.869
## MMPI_HY   |    1.522 |   1.522 |    1 |     2.422 |   0.120 | 0.000 |         0.001 |   0.000 |           0.000 |     0.000 |    0.023 | 0.343
## MMPI_PD   |   38.312 |  38.312 |    1 |    60.946 |  < .001 | 0.011 |         0.014 |   0.011 |           0.013 |     0.011 |    0.117 | 1.000
## MMPI_MF   |   53.998 |  53.998 |    1 |    85.900 |  < .001 | 0.015 |         0.019 |   0.015 |           0.019 |     0.015 |    0.139 | 1.000
## MMPI_PA   |    5.786 |   5.786 |    1 |     9.205 |   0.002 | 0.002 |         0.002 |   0.001 |           0.002 |     0.001 |    0.046 | 0.859
## MMPI_PT   |    2.958 |   2.958 |    1 |     4.705 |   0.030 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.033 | 0.583
## MMPI_SC   |   10.214 |  10.214 |    1 |    16.248 |  < .001 | 0.003 |         0.004 |   0.003 |           0.003 |     0.003 |    0.060 | 0.981
## MMPI_MA   |    8.259 |   8.259 |    1 |    13.138 |  < .001 | 0.002 |         0.003 |   0.002 |           0.003 |     0.002 |    0.054 | 0.952
## MMPI_SI   |   13.249 |  13.249 |    1 |    21.077 |  < .001 | 0.004 |         0.005 |   0.004 |           0.004 |     0.004 |    0.069 | 0.996
## MMPI_ES   |    1.283 |   1.283 |    1 |     2.041 |   0.153 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.021 | 0.298
## Residuals | 2793.549 |   0.629 | 4444 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 2 × 3
##   group sum_omegasq sum_omega
##   <chr>       <dbl>     <dbl>
## 1 g           0.16      0.4  
## 2 MMPI        0.047     0.217

compare_total_effects(
  lm(as.formula(str_glue("unemployment ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), data = d),
  groups1
)

## $anova_stats
## term      |    sumsq | meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ---------------------------------------------------------------------------------------------------------------------------------------------
## g         |   43.907 | 43.907 |    1 |    41.322 |  < .001 | 0.009 |         0.009 |   0.009 |           0.009 |     0.009 |    0.097 | 1.000
## MMPI_L    |    8.392 |  8.392 |    1 |     7.898 |   0.005 | 0.002 |         0.002 |   0.002 |           0.002 |     0.002 |    0.042 | 0.802
## MMPI_F    |   19.409 | 19.409 |    1 |    18.266 |  < .001 | 0.004 |         0.004 |   0.004 |           0.004 |     0.004 |    0.064 | 0.990
## MMPI_K    |    4.193 |  4.193 |    1 |     3.946 |   0.047 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.030 | 0.511
## MMPI_HS   |    0.773 |  0.773 |    1 |     0.727 |   0.394 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.013 | 0.137
## MMPI_D    |   12.217 | 12.217 |    1 |    11.498 |   0.001 | 0.003 |         0.003 |   0.002 |           0.002 |     0.002 |    0.051 | 0.924
## MMPI_HY   |    0.448 |  0.448 |    1 |     0.422 |   0.516 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.010 | 0.100
## MMPI_PD   |   33.274 | 33.274 |    1 |    31.315 |  < .001 | 0.007 |         0.007 |   0.007 |           0.007 |     0.007 |    0.084 | 1.000
## MMPI_MF   |    0.099 |  0.099 |    1 |     0.093 |   0.761 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.005 | 0.061
## MMPI_PA   |    0.298 |  0.298 |    1 |     0.280 |   0.597 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.008 | 0.083
## MMPI_PT   |   14.259 | 14.259 |    1 |    13.420 |  < .001 | 0.003 |         0.003 |   0.003 |           0.003 |     0.003 |    0.055 | 0.956
## MMPI_SC   |    8.188 |  8.188 |    1 |     7.706 |   0.006 | 0.002 |         0.002 |   0.001 |           0.002 |     0.001 |    0.042 | 0.793
## MMPI_MA   |    1.707 |  1.707 |    1 |     1.607 |   0.205 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.019 | 0.245
## MMPI_SI   |   10.342 | 10.342 |    1 |     9.734 |   0.002 | 0.002 |         0.002 |   0.002 |           0.002 |     0.002 |    0.047 | 0.877
## MMPI_ES   |    7.930 |  7.930 |    1 |     7.463 |   0.006 | 0.002 |         0.002 |   0.001 |           0.001 |     0.001 |    0.041 | 0.780
## Residuals | 4708.191 |  1.063 | 4431 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 2 × 3
##   group sum_omegasq sum_omega
##   <chr>       <dbl>     <dbl>
## 1 g           0.009    0.0949
## 2 MMPI        0.023    0.152

compare_total_effects(
  lm(as.formula(str_glue("military_rank ~ g + {str_c(MMPI_vars, collapse = ' + ')}")), data = d),
  groups1
)

## $anova_stats
## term      |    sumsq |  meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ----------------------------------------------------------------------------------------------------------------------------------------------
## g         |  159.978 | 159.978 |    1 |   155.954 |  < .001 | 0.033 |         0.034 |   0.033 |           0.034 |     0.033 |    0.187 | 1.000
## MMPI_L    |    9.115 |   9.115 |    1 |     8.886 |   0.003 | 0.002 |         0.002 |   0.002 |           0.002 |     0.002 |    0.045 | 0.846
## MMPI_F    |    8.319 |   8.319 |    1 |     8.110 |   0.004 | 0.002 |         0.002 |   0.002 |           0.002 |     0.002 |    0.043 | 0.813
## MMPI_K    |    8.287 |   8.287 |    1 |     8.078 |   0.005 | 0.002 |         0.002 |   0.001 |           0.002 |     0.001 |    0.043 | 0.811
## MMPI_HS   |    1.059 |   1.059 |    1 |     1.032 |   0.310 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.015 | 0.174
## MMPI_D    |    0.107 |   0.107 |    1 |     0.105 |   0.746 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.005 | 0.062
## MMPI_HY   |    9.021 |   9.021 |    1 |     8.794 |   0.003 | 0.002 |         0.002 |   0.002 |           0.002 |     0.002 |    0.044 | 0.843
## MMPI_PD   |   79.707 |  79.707 |    1 |    77.702 |  < .001 | 0.016 |         0.017 |   0.016 |           0.017 |     0.016 |    0.132 | 1.000
## MMPI_MF   |    2.147 |   2.147 |    1 |     2.093 |   0.148 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.022 | 0.304
## MMPI_PA   |    2.869 |   2.869 |    1 |     2.796 |   0.095 | 0.001 |         0.001 |   0.000 |           0.000 |     0.000 |    0.025 | 0.387
## MMPI_PT   |    0.672 |   0.672 |    1 |     0.655 |   0.418 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.012 | 0.128
## MMPI_SC   |    0.002 |   0.002 |    1 |     0.002 |   0.967 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.001 | 0.050
## MMPI_MA   |    1.680 |   1.680 |    1 |     1.638 |   0.201 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.019 | 0.249
## MMPI_SI   |    2.614 |   2.614 |    1 |     2.548 |   0.110 | 0.001 |         0.001 |   0.000 |           0.000 |     0.000 |    0.024 | 0.358
## MMPI_ES   |    0.077 |   0.077 |    1 |     0.075 |   0.784 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.004 | 0.059
## Residuals | 4560.725 |   1.026 | 4446 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 2 × 3
##   group sum_omegasq sum_omega
##   <chr>       <dbl>     <dbl>
## 1 g           0.033     0.182
## 2 MMPI        0.023     0.152

#control for race and age
groups2 = list(g = "g", MMPI = MMPI_vars, age_race = c("age", "race"))
compare_total_effects(
  lm(as.formula(str_glue("income ~ g + {str_c(MMPI_vars, collapse = ' + ')} + age + race")), data = d),
  groups2
)

## $anova_stats
## term      |    sumsq |  meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ----------------------------------------------------------------------------------------------------------------------------------------------
## g         |  210.325 | 210.325 |    1 |   256.456 |  < .001 | 0.053 |         0.056 |   0.053 |           0.055 |     0.053 |    0.243 | 1.000
## MMPI_L    |   13.489 |  13.489 |    1 |    16.448 |  < .001 | 0.003 |         0.004 |   0.003 |           0.004 |     0.003 |    0.061 | 0.982
## MMPI_F    |    2.039 |   2.039 |    1 |     2.486 |   0.115 | 0.001 |         0.001 |   0.000 |           0.000 |     0.000 |    0.024 | 0.351
## MMPI_K    |    2.738 |   2.738 |    1 |     3.339 |   0.068 | 0.001 |         0.001 |   0.000 |           0.001 |     0.000 |    0.028 | 0.447
## MMPI_HS   |    1.380 |   1.380 |    1 |     1.683 |   0.195 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.020 | 0.254
## MMPI_D    |    0.399 |   0.399 |    1 |     0.486 |   0.486 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.011 | 0.107
## MMPI_HY   |    4.336 |   4.336 |    1 |     5.287 |   0.022 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.035 | 0.633
## MMPI_PD   |   43.811 |  43.811 |    1 |    53.420 |  < .001 | 0.011 |         0.012 |   0.011 |           0.012 |     0.011 |    0.111 | 1.000
## MMPI_MF   |    1.123 |   1.123 |    1 |     1.370 |   0.242 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.018 | 0.216
## MMPI_PA   |    1.249 |   1.249 |    1 |     1.523 |   0.217 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.019 | 0.235
## MMPI_PT   |    4.802 |   4.802 |    1 |     5.855 |   0.016 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.037 | 0.677
## MMPI_SC   |    1.256 |   1.256 |    1 |     1.532 |   0.216 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.019 | 0.236
## MMPI_MA   |    5.066 |   5.066 |    1 |     6.178 |   0.013 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.038 | 0.700
## MMPI_SI   |    0.026 |   0.026 |    1 |     0.032 |   0.859 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.003 | 0.054
## MMPI_ES   |   16.996 |  16.996 |    1 |    20.724 |  < .001 | 0.004 |         0.005 |   0.004 |           0.004 |     0.004 |    0.069 | 0.995
## age       |   73.667 |  73.667 |    1 |    89.825 |  < .001 | 0.019 |         0.020 |   0.018 |           0.020 |     0.018 |    0.144 | 1.000
## race      |    6.621 |   1.655 |    4 |     2.018 |   0.089 | 0.002 |         0.002 |   0.001 |           0.001 |     0.001 |    0.043 | 0.609
## Residuals | 3571.637 |   0.820 | 4355 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 3 × 3
##   group    sum_omegasq sum_omega
##   <chr>          <dbl>     <dbl>
## 1 g              0.053     0.230
## 2 MMPI           0.021     0.145
## 3 age_race       0.019     0.138

compare_total_effects(
  lm(as.formula(str_glue("education ~ g + {str_c(MMPI_vars, collapse = ' + ')} + age + race")), data = d),
  groups2
)

## $anova_stats
## term      |    sumsq |  meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ----------------------------------------------------------------------------------------------------------------------------------------------
## g         |  638.306 | 638.306 |    1 |  1072.575 |  < .001 | 0.180 |         0.195 |   0.180 |           0.194 |     0.180 |    0.492 | 1.000
## MMPI_L    |    1.001 |   1.001 |    1 |     1.681 |   0.195 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.019 | 0.254
## MMPI_F    |    1.422 |   1.422 |    1 |     2.389 |   0.122 | 0.000 |         0.001 |   0.000 |           0.000 |     0.000 |    0.023 | 0.340
## MMPI_K    |   12.638 |  12.638 |    1 |    21.236 |  < .001 | 0.004 |         0.005 |   0.003 |           0.005 |     0.003 |    0.069 | 0.996
## MMPI_HS   |    4.197 |   4.197 |    1 |     7.053 |   0.008 | 0.001 |         0.002 |   0.001 |           0.001 |     0.001 |    0.040 | 0.757
## MMPI_D    |    2.705 |   2.705 |    1 |     4.546 |   0.033 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.032 | 0.568
## MMPI_HY   |    0.250 |   0.250 |    1 |     0.420 |   0.517 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.010 | 0.099
## MMPI_PD   |   28.413 |  28.413 |    1 |    47.744 |  < .001 | 0.008 |         0.011 |   0.008 |           0.010 |     0.008 |    0.104 | 1.000
## MMPI_MF   |   42.522 |  42.522 |    1 |    71.451 |  < .001 | 0.012 |         0.016 |   0.012 |           0.016 |     0.012 |    0.127 | 1.000
## MMPI_PA   |    5.876 |   5.876 |    1 |     9.874 |   0.002 | 0.002 |         0.002 |   0.001 |           0.002 |     0.001 |    0.047 | 0.881
## MMPI_PT   |    0.429 |   0.429 |    1 |     0.721 |   0.396 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.013 | 0.136
## MMPI_SC   |    5.958 |   5.958 |    1 |    10.012 |   0.002 | 0.002 |         0.002 |   0.002 |           0.002 |     0.002 |    0.047 | 0.886
## MMPI_MA   |    3.301 |   3.301 |    1 |     5.547 |   0.019 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.035 | 0.654
## MMPI_SI   |   10.856 |  10.856 |    1 |    18.242 |  < .001 | 0.003 |         0.004 |   0.003 |           0.004 |     0.003 |    0.064 | 0.990
## MMPI_ES   |    0.001 |   0.001 |    1 |     0.001 |   0.973 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.001 | 0.050
## age       |   48.973 |  48.973 |    1 |    82.292 |  < .001 | 0.014 |         0.018 |   0.014 |           0.018 |     0.014 |    0.136 | 1.000
## race      |   96.416 |  24.104 |    4 |    40.503 |  < .001 | 0.027 |         0.035 |   0.027 |           0.034 |     0.027 |    0.191 | 1.000
## Residuals | 2641.717 |   0.595 | 4439 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 3 × 3
##   group    sum_omegasq sum_omega
##   <chr>          <dbl>     <dbl>
## 1 g              0.18      0.424
## 2 MMPI           0.032     0.179
## 3 age_race       0.041     0.202

compare_total_effects(
  lm(as.formula(str_glue("unemployment ~ g + {str_c(MMPI_vars, collapse = ' + ')} + age + race")), data = d),
  groups2
)

## $anova_stats
## term      |    sumsq | meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ---------------------------------------------------------------------------------------------------------------------------------------------
## g         |   20.962 | 20.962 |    1 |    19.943 |  < .001 | 0.004 |         0.004 |   0.004 |           0.004 |     0.004 |    0.067 | 0.994
## MMPI_L    |    7.377 |  7.377 |    1 |     7.018 |   0.008 | 0.002 |         0.002 |   0.001 |           0.001 |     0.001 |    0.040 | 0.755
## MMPI_F    |   21.525 | 21.525 |    1 |    20.479 |  < .001 | 0.004 |         0.005 |   0.004 |           0.004 |     0.004 |    0.068 | 0.995
## MMPI_K    |    4.074 |  4.074 |    1 |     3.876 |   0.049 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.030 | 0.504
## MMPI_HS   |    0.688 |  0.688 |    1 |     0.654 |   0.419 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.012 | 0.128
## MMPI_D    |   12.279 | 12.279 |    1 |    11.682 |   0.001 | 0.003 |         0.003 |   0.002 |           0.002 |     0.002 |    0.051 | 0.928
## MMPI_HY   |    1.598 |  1.598 |    1 |     1.520 |   0.218 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.019 | 0.234
## MMPI_PD   |   25.135 | 25.135 |    1 |    23.914 |  < .001 | 0.005 |         0.005 |   0.005 |           0.005 |     0.005 |    0.074 | 0.998
## MMPI_MF   |    0.722 |  0.722 |    1 |     0.686 |   0.407 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.012 | 0.132
## MMPI_PA   |    0.028 |  0.028 |    1 |     0.026 |   0.871 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.002 | 0.053
## MMPI_PT   |   13.309 | 13.309 |    1 |    12.662 |  < .001 | 0.003 |         0.003 |   0.003 |           0.003 |     0.003 |    0.053 | 0.945
## MMPI_SC   |    6.925 |  6.925 |    1 |     6.588 |   0.010 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.039 | 0.728
## MMPI_MA   |    0.311 |  0.311 |    1 |     0.296 |   0.586 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.008 | 0.085
## MMPI_SI   |    8.552 |  8.552 |    1 |     8.137 |   0.004 | 0.002 |         0.002 |   0.002 |           0.002 |     0.002 |    0.043 | 0.814
## MMPI_ES   |    6.460 |  6.460 |    1 |     6.146 |   0.013 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.037 | 0.698
## age       |   18.943 | 18.943 |    1 |    18.023 |  < .001 | 0.004 |         0.004 |   0.004 |           0.004 |     0.004 |    0.064 | 0.989
## race      |   38.327 |  9.582 |    4 |     9.116 |  < .001 | 0.008 |         0.008 |   0.007 |           0.007 |     0.007 |    0.091 | 1.000
## Residuals | 4652.076 |  1.051 | 4426 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 3 × 3
##   group    sum_omegasq sum_omega
##   <chr>          <dbl>     <dbl>
## 1 g              0.004    0.0632
## 2 MMPI           0.02     0.141 
## 3 age_race       0.011    0.105

compare_total_effects(
  lm(as.formula(str_glue("military_rank ~ g + {str_c(MMPI_vars, collapse = ' + ')} + age + race")), data = d),
  groups2
)

## $anova_stats
## term      |    sumsq |  meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ----------------------------------------------------------------------------------------------------------------------------------------------
## g         |  101.845 | 101.845 |    1 |   105.113 |  < .001 | 0.021 |         0.023 |   0.021 |           0.023 |     0.021 |    0.154 | 1.000
## MMPI_L    |    5.489 |   5.489 |    1 |     5.665 |   0.017 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.036 | 0.663
## MMPI_F    |    7.621 |   7.621 |    1 |     7.866 |   0.005 | 0.002 |         0.002 |   0.001 |           0.002 |     0.001 |    0.042 | 0.801
## MMPI_K    |    4.647 |   4.647 |    1 |     4.796 |   0.029 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.033 | 0.591
## MMPI_HS   |    0.812 |   0.812 |    1 |     0.838 |   0.360 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.014 | 0.150
## MMPI_D    |    0.020 |   0.020 |    1 |     0.021 |   0.885 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.002 | 0.052
## MMPI_HY   |    4.025 |   4.025 |    1 |     4.154 |   0.042 | 0.001 |         0.001 |   0.001 |           0.001 |     0.001 |    0.031 | 0.531
## MMPI_PD   |   51.963 |  51.963 |    1 |    53.631 |  < .001 | 0.011 |         0.012 |   0.011 |           0.012 |     0.011 |    0.110 | 1.000
## MMPI_MF   |    1.120 |   1.120 |    1 |     1.156 |   0.282 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.016 | 0.189
## MMPI_PA   |    1.061 |   1.061 |    1 |     1.095 |   0.295 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.016 | 0.182
## MMPI_PT   |    1.287 |   1.287 |    1 |     1.328 |   0.249 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.017 | 0.211
## MMPI_SC   |    0.070 |   0.070 |    1 |     0.073 |   0.787 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.004 | 0.058
## MMPI_MA   |    0.106 |   0.106 |    1 |     0.109 |   0.741 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.005 | 0.063
## MMPI_SI   |    0.939 |   0.939 |    1 |     0.969 |   0.325 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.015 | 0.166
## MMPI_ES   |    0.072 |   0.072 |    1 |     0.075 |   0.785 | 0.000 |         0.000 |   0.000 |           0.000 |     0.000 |    0.004 | 0.059
## age       |  216.309 | 216.309 |    1 |   223.250 |  < .001 | 0.046 |         0.048 |   0.045 |           0.047 |     0.045 |    0.224 | 1.000
## race      |   44.513 |  11.128 |    4 |    11.485 |  < .001 | 0.009 |         0.010 |   0.009 |           0.009 |     0.009 |    0.102 | 1.000
## Residuals | 4302.915 |   0.969 | 4441 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 3 × 3
##   group    sum_omegasq sum_omega
##   <chr>          <dbl>     <dbl>
## 1 g              0.021     0.145
## 2 MMPI           0.015     0.122
## 3 age_race       0.054     0.232

MMPI g vs. g

#predict income
list(
  ols(income ~ g, data = d),
  ols(income ~ MMPI_p, data = d),
  ols(income ~ g + MMPI_p, data = d)
  ) %>% 
  summarize_models()

#education
list(
  ols(education ~ g, data = d),
  ols(education ~ MMPI_p, data = d),
  ols(education ~ g + MMPI_p, data = d)
  ) %>% 
  summarize_models()

#occupational status
# list(
#   ols(occu_status ~ g, data = d),
#   ols(occu_status ~ MMPI_p, data = d),
#   ols(occu_status ~ g + MMPI_p, data = d)
#   ) %>% 
#   summarize_models()

#unemployment
list(
  ols(unemployment ~ g, data = d),
  ols(unemployment ~ MMPI_p, data = d),
  ols(unemployment ~ g + MMPI_p, data = d)
  ) %>% 
  summarize_models()

#military rank
list(
  ols(military_rank ~ g, data = d),
  ols(military_rank ~ MMPI_p, data = d),
  ols(military_rank ~ g + MMPI_p, data = d)
  ) %>% 
  summarize_models()

#discharge status
list(
  lrm(discharge_ok ~ g, data = d),
  lrm(discharge_ok ~ MMPI_p, data = d),
  lrm(discharge_ok ~ g + MMPI_p, data = d)
  )

## [[1]]
## Logistic Regression Model
##  
##  lrm(formula = discharge_ok ~ g, data = d)
##  
##                         Model Likelihood      Discrimination    Rank Discrim.    
##                               Ratio Test             Indexes          Indexes    
##  Obs          4462    LR chi2      39.50      R2       0.052    C       0.701    
##   FALSE         84    d.f.             1     R2(1,4462)0.009    Dxy     0.402    
##   TRUE        4378    Pr(> chi2) <0.0001    R2(1,247.3)0.144    gamma   0.404    
##  max |deriv| 2e-14                            Brier    0.018    tau-a   0.015    
##  
##            Coef   S.E.   Wald Z Pr(>|Z|)
##  Intercept 4.3110 0.1461 29.51  <0.0001 
##  g         0.6463 0.1051  6.15  <0.0001 
##  
## 
## [[2]]
## Logistic Regression Model
##  
##  lrm(formula = discharge_ok ~ MMPI_p, data = d)
##  
##                         Model Likelihood      Discrimination    Rank Discrim.    
##                               Ratio Test             Indexes          Indexes    
##  Obs          4462    LR chi2      16.80      R2       0.022    C       0.642    
##   FALSE         84    d.f.             1     R2(1,4462)0.004    Dxy     0.284    
##   TRUE        4378    Pr(> chi2) <0.0001    R2(1,247.3)0.062    gamma   0.286    
##  max |deriv| 4e-09                            Brier    0.018    tau-a   0.010    
##  
##            Coef    S.E.   Wald Z Pr(>|Z|)
##  Intercept  4.0841 0.1238 32.99  <0.0001 
##  MMPI_p    -0.4552 0.1119 -4.07  <0.0001 
##  
## 
## [[3]]
## Logistic Regression Model
##  
##  lrm(formula = discharge_ok ~ g + MMPI_p, data = d)
##  
##                         Model Likelihood      Discrimination    Rank Discrim.    
##                               Ratio Test             Indexes          Indexes    
##  Obs          4462    LR chi2      43.48      R2       0.057    C       0.712    
##   FALSE         84    d.f.             2     R2(2,4462)0.009    Dxy     0.424    
##   TRUE        4378    Pr(> chi2) <0.0001    R2(2,247.3)0.154    gamma   0.426    
##  max |deriv| 2e-13                            Brier    0.018    tau-a   0.016    
##  
##            Coef    S.E.   Wald Z Pr(>|Z|)
##  Intercept  4.3362 0.1482 29.25  <0.0001 
##  g          0.5682 0.1120  5.07  <0.0001 
##  MMPI_p    -0.2384 0.1200 -1.99  0.0469  
##

#marrital status
list(
  lrm(married ~ g, data = d),
  lrm(married ~ MMPI_p, data = d),
  lrm(married ~ g + MMPI_p, data = d)
  )

## [[1]]
## Frequencies of Missing Values Due to Each Variable
## married       g 
##       3       0 
## 
## Logistic Regression Model
##  
##  lrm(formula = married ~ g, data = d)
##  
##  
##                         Model Likelihood       Discrimination    Rank Discrim.    
##                               Ratio Test              Indexes          Indexes    
##  Obs          4459    LR chi2      44.35       R2       0.014    C       0.561    
##   FALSE       1184    d.f.             1      R2(1,4459)0.010    Dxy     0.123    
##   TRUE        3275    Pr(> chi2) <0.0001    R2(1,2608.8)0.016    gamma   0.123    
##  max |deriv| 4e-14                             Brier    0.193    tau-a   0.048    
##  
##            Coef   S.E.   Wald Z Pr(>|Z|)
##  Intercept 1.0698 0.0354 30.20  <0.0001 
##  g         0.2107 0.0318  6.64  <0.0001 
##  
## 
## [[2]]
## Frequencies of Missing Values Due to Each Variable
## married  MMPI_p 
##       3       0 
## 
## Logistic Regression Model
##  
##  lrm(formula = married ~ MMPI_p, data = d)
##  
##  
##                         Model Likelihood       Discrimination    Rank Discrim.    
##                               Ratio Test              Indexes          Indexes    
##  Obs          4459    LR chi2      95.30       R2       0.031    C       0.597    
##   FALSE       1184    d.f.             1      R2(1,4459)0.021    Dxy     0.193    
##   TRUE        3275    Pr(> chi2) <0.0001    R2(1,2608.8)0.036    gamma   0.193    
##  max |deriv| 2e-13                             Brier    0.191    tau-a   0.075    
##  
##            Coef    S.E.   Wald Z Pr(>|Z|)
##  Intercept  1.0667 0.0351 30.35  <0.0001 
##  MMPI_p    -0.3361 0.0350 -9.59  <0.0001 
##  
## 
## [[3]]
## Frequencies of Missing Values Due to Each Variable
## married       g  MMPI_p 
##       3       0       0 
## 
## Logistic Regression Model
##  
##  lrm(formula = married ~ g + MMPI_p, data = d)
##  
##  
##                         Model Likelihood       Discrimination    Rank Discrim.    
##                               Ratio Test              Indexes          Indexes    
##  Obs          4459    LR chi2     107.37       R2       0.035    C       0.602    
##   FALSE       1184    d.f.             2      R2(2,4459)0.023    Dxy     0.204    
##   TRUE        3275    Pr(> chi2) <0.0001    R2(2,2608.8)0.040    gamma   0.204    
##  max |deriv| 3e-13                             Brier    0.190    tau-a   0.080    
##  
##            Coef    S.E.   Wald Z Pr(>|Z|)
##  Intercept  1.0894 0.0360 30.28  <0.0001 
##  g          0.1178 0.0339  3.47  0.0005  
##  MMPI_p    -0.2919 0.0372 -7.84  <0.0001 
##

#control for race and age
groups3 = list(g = "g", MMPI = "MMPI_p", age_race = c("age", "race"))
compare_total_effects(
  lm(as.formula(str_glue("income ~ g + MMPI_p + age + race")), data = d),
  groups3
)

## $anova_stats
## term      |    sumsq |  meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ----------------------------------------------------------------------------------------------------------------------------------------------
## g         |  354.070 | 354.070 |    1 |   420.604 |  < .001 | 0.084 |         0.088 |   0.084 |           0.087 |     0.084 |    0.310 | 1.000
## MMPI_p    |   95.433 |  95.433 |    1 |   113.367 |  < .001 | 0.023 |         0.025 |   0.022 |           0.025 |     0.022 |    0.161 | 1.000
## age       |   83.352 |  83.352 |    1 |    99.015 |  < .001 | 0.020 |         0.022 |   0.020 |           0.022 |     0.020 |    0.151 | 1.000
## race      |    9.299 |   2.325 |    4 |     2.762 |   0.026 | 0.002 |         0.003 |   0.001 |           0.002 |     0.001 |    0.050 | 0.764
## Residuals | 3677.038 |   0.842 | 4368 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 3 × 3
##   group    sum_omegasq sum_omega
##   <chr>          <dbl>     <dbl>
## 1 g              0.084     0.290
## 2 MMPI           0.022     0.148
## 3 age_race       0.021     0.145

compare_total_effects(
  lm(as.formula(str_glue("education ~ g + MMPI_p + age + race")), data = d),
  groups3
)

## $anova_stats
## term      |    sumsq |   meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## -----------------------------------------------------------------------------------------------------------------------------------------------
## g         | 1102.071 | 1102.071 |    1 |  1772.141 |  < .001 | 0.269 |         0.285 |   0.269 |           0.284 |     0.269 |    0.631 |     1
## MMPI_p    |   27.874 |   27.874 |    1 |    44.822 |  < .001 | 0.007 |         0.010 |   0.007 |           0.010 |     0.007 |    0.100 |     1
## age       |   49.151 |   49.151 |    1 |    79.036 |  < .001 | 0.012 |         0.017 |   0.012 |           0.017 |     0.012 |    0.133 |     1
## race      |  148.644 |   37.161 |    4 |    59.755 |  < .001 | 0.036 |         0.051 |   0.036 |           0.050 |     0.036 |    0.232 |     1
## Residuals | 2768.639 |    0.622 | 4452 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 3 × 3
##   group    sum_omegasq sum_omega
##   <chr>          <dbl>     <dbl>
## 1 g              0.269    0.519 
## 2 MMPI           0.007    0.0837
## 3 age_race       0.048    0.219

compare_total_effects(
  lm(as.formula(str_glue("unemployment ~ g + MMPI_p + age + race")), data = d),
  groups3
)

## $anova_stats
## term      |    sumsq |  meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ----------------------------------------------------------------------------------------------------------------------------------------------
## g         |   53.008 |  53.008 |    1 |    47.428 |  < .001 | 0.010 |         0.011 |   0.010 |           0.010 |     0.010 |    0.103 | 1.000
## MMPI_p    |  160.674 | 160.674 |    1 |   143.759 |  < .001 | 0.031 |         0.031 |   0.030 |           0.031 |     0.030 |    0.180 | 1.000
## age       |   30.577 |  30.577 |    1 |    27.358 |  < .001 | 0.006 |         0.006 |   0.006 |           0.006 |     0.006 |    0.079 | 0.999
## race      |   48.542 |  12.136 |    4 |    10.858 |  < .001 | 0.009 |         0.010 |   0.008 |           0.009 |     0.008 |    0.099 | 1.000
## Residuals | 4961.305 |   1.118 | 4439 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 3 × 3
##   group    sum_omegasq sum_omega
##   <chr>          <dbl>     <dbl>
## 1 g              0.01      0.1  
## 2 MMPI           0.03      0.173
## 3 age_race       0.014     0.118

compare_total_effects(
  lm(as.formula(str_glue("military_rank ~ g + MMPI_p + age + race")), data = d),
  groups3
)

## $anova_stats
## term      |    sumsq |  meansq |   df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
## ----------------------------------------------------------------------------------------------------------------------------------------------
## g         |  145.557 | 145.557 |    1 |   147.437 |  < .001 | 0.030 |         0.032 |   0.030 |           0.032 |     0.030 |    0.182 |     1
## MMPI_p    |   31.670 |  31.670 |    1 |    32.079 |  < .001 | 0.006 |         0.007 |   0.006 |           0.007 |     0.006 |    0.085 |     1
## age       |  260.625 | 260.625 |    1 |   263.991 |  < .001 | 0.053 |         0.056 |   0.053 |           0.056 |     0.053 |    0.243 |     1
## race      |   57.822 |  14.456 |    4 |    14.642 |  < .001 | 0.012 |         0.013 |   0.011 |           0.012 |     0.011 |    0.115 |     1
## Residuals | 4397.205 |   0.987 | 4454 |           |         |       |               |         |                 |           |          |      
## 
## $sums
## # A tibble: 3 × 3
##   group    sum_omegasq sum_omega
##   <chr>          <dbl>     <dbl>
## 1 g              0.03     0.173 
## 2 MMPI           0.006    0.0775
## 3 age_race       0.064    0.253

VES: personality vs. g for success

22 October, 2022

Init

Functions

Data

Recode

Analysis

MMPI scales vs. g

MMPI g vs. g