quick-introduction

Version Note: Up-to-date with v0.5.0

library(psycModel)

Why would you want to use this package?

TLDR:
1) It’s a beginner-friendly R package for statistical analysis in social science.
2) Fitting models, plotting, checking the goodness of fit, and identifying model assumption violations all in one place.
3) Beautiful and easy-to-read output. Check out this example now.

Some Examples:

Model Summary

The model_summary function will produce all of the relevant test statistics for regression models. See an example below.

mod_1 = lm(data = iris, Petal.Length ~ Petal.Width*Sepal.Length)
# you can also check assumption_plot by setting assumption_plot = TRUE
model_summary(mod_1,assumption_plot = F) 

 
Model Summary
Model Type = Linear regression
Outcome = Petal.Length
Predictors = Petal.Width, Sepal.Length

Model Estimates
────────────────────────────────────────────────────────────────────────────────────────
                 Parameter  Coefficient     SE       t   df          p            95% CI
────────────────────────────────────────────────────────────────────────────────────────
               (Intercept)       -3.248  0.596  -5.451  146  0.000 ***  [-4.426, -2.070]
               Petal.Width        2.971  0.358   8.291  146  0.000 ***  [ 2.263,  3.679]
              Sepal.Length        0.876  0.117   7.504  146  0.000 ***  [ 0.645,  1.106]
  Petal.Width:Sepal.Length       -0.222  0.064  -3.485  146  0.001 ***  [-0.349, -0.096]
────────────────────────────────────────────────────────────────────────────────────────
*** p < 0.001, ** p < 0.01, * p < 0.05, + p < 0.1

Goodness of Fit
─────────────────────────────────────────────────────────────
      AIC     AICc      BIC     R²  R²_adjusted   RMSE      σ
─────────────────────────────────────────────────────────────
  148.192  148.608  163.245  0.952        0.952  0.384  0.389
─────────────────────────────────────────────────────────────

Interaction Plot

Since our models have an interaction, we may want to visualize it. Let’s create an interaction plot first. You can modify the plot using some of the settings. You can also try out the polynomial regression plots and the ANOVA plots.

interaction_plot(mod_1,verbose = F) # verbose is set to TRUE by default to show the plot data. 

# You can also modify the 
interaction_plot(mod_1,
                 response_var_name = 'Petal Length', # you can rename the variable names (e.g., getting rid of the dot)
                 predict_var1_name = 'Petal Width',
                 predict_var2_name = 'Sepal Length',
                 predict_var1_level = c(0.44,1.19,1.96), # you may want to show the mean of the variable
                 predict_var1_level_name = c('-1 SD','Mean','+1 SD'), # you may also want to adjust the labels to be more intuitive
                 verbose = F)

Simple Slopes

After obtaining the interaction plot, you may also want to get the simple slopes of the interaction.

simple_slope(mod_1)
$simple_slope_df
  Sepal.Length Level     Est.       S.E. ci.lower ci.upper   t val.
1                Low 1.855359 0.07887932 1.699467 2.011252 23.52149
2               Mean 1.671132 0.07591186 1.521104 1.821160 22.01411
3               High 1.486904 0.10436269 1.280648 1.693161 14.24747
             p
1 1.615129e-51
2 3.090242e-48
3 2.027575e-29

$jn_plot

Descriptive Table

This package can also help you in preparing a table that includes means, standard deviations, and correlations. For additional options, refer to ?descriptive_table.

descriptive_table(iris,c(Petal.Width,Sepal.Length,Petal.Length))
Model Summary
Model Type = Correlation
Model Method = pearson
Adjustment Method = none

─────────────────────────────────────────
           Var  Petal.Width  Sepal.Length
─────────────────────────────────────────
   Petal.Width                           
  Sepal.Length    0.818 ***              
  Petal.Length    0.963 ***     0.872 ***
─────────────────────────────────────────
Note: * p < 0.05, ** p < 0.01, *** p < 0.001

Model Summary
Model Type = Descriptive Statistics

───────────────────────────────────────────────────────
           Var   mean     sd  Petal.Width  Sepal.Length
───────────────────────────────────────────────────────
   Petal.Width  1.199  0.762                           
  Sepal.Length  5.843  0.828    0.818 ***              
  Petal.Length  3.758  1.765    0.963 ***     0.872 ***
───────────────────────────────────────────────────────
descriptive_table(iris,c(Petal.Width,Sepal.Length,Petal.Length),descriptive_indicator = c('mean','sd','cor','missing','kurtosis')) # you can request more parameters optionally
Model Summary
Model Type = Correlation
Model Method = pearson
Adjustment Method = none

─────────────────────────────────────────
           Var  Petal.Width  Sepal.Length
─────────────────────────────────────────
   Petal.Width                           
  Sepal.Length    0.818 ***              
  Petal.Length    0.963 ***     0.872 ***
─────────────────────────────────────────
Note: * p < 0.05, ** p < 0.01, *** p < 0.001

Model Summary
Model Type = Descriptive Statistics

────────────────────────────────────────────────────────────────────────────
           Var  missing_n   mean     sd  kurtosis  Petal.Width  Sepal.Length
────────────────────────────────────────────────────────────────────────────
   Petal.Width          0  1.199  0.762    -1.358                           
  Sepal.Length          0  5.843  0.828    -0.606    0.818 ***              
  Petal.Length          0  3.758  1.765    -1.417    0.963 ***     0.872 ***
────────────────────────────────────────────────────────────────────────────

Cronbach alpha

You can get the Cronbach alphas very simply (it will call the psych::alpha() function). If you need, you can also get separate Cronbach alphas by groups (e.g., when using multilevel analyses).

cronbach_alpha(x1:x3,x4:x6,x7:x9, 
               var_name = c('visual','textual','verbal'),
               data = lavaan::HolzingerSwineford1939)

 
Model Summary
Model Type = Cronbach Alpha Reliability Analysis
Model Specification: 
 visual = x1 + x2 + x3 
 textual = x4 + x5 + x6 
 verbal = x7 + x8 + x9 

───────────────────────────────
      Var  raw_alpha  std_alpha
───────────────────────────────
   visual      0.626      0.627
  textual      0.883      0.885
   verbal      0.688      0.690
───────────────────────────────
cronbach_alpha(x1:x3,x4:x6,x7:x9, 
               var_name = c('visual','textual','verbal'),
               group = 'sex',
               data = lavaan::HolzingerSwineford1939)

 
Model Summary
Model Type = Cronbach Alpha Reliability Analysis
Model Specification: 
 visual = x1 + x2 + x3 
 textual = x4 + x5 + x6 
 verbal = x7 + x8 + x9 

────────────────────────────────────
      Var  sex  raw_alpha  std_alpha
────────────────────────────────────
   visual    1      0.568      0.572
   visual    2      0.664      0.663
  textual    1      0.872      0.874
  textual    2      0.892      0.895
   verbal    1      0.697      0.693
   verbal    2      0.686      0.697
────────────────────────────────────

Confirmatory Factor Analysis

CFA model is fitted using lavaan::cfa(). You can pass multiple factor (in the below example, x1, x2, x3 represent one factor, x4,x5,x6 represent another factor etc.). It will show you the fit measure, factor loading, and goodness of fit based on cut-off criteria (you should review literature for the cut-off criteria as the recommendations are subjected to changes). You can also try measurement_invariance().

cfa_summary(
   data = lavaan::HolzingerSwineford1939,
   x1:x3,
   x4:x6,
   x7:x9
 )

 
Model Summary
Model Type = Confirmatory Factor Analysis
Estimator: ML
Model Formula = 
. DV1 =~ x1 + x2 + x3
  DV2 =~ x4 + x5 + x6
  DV3 =~ x7 + x8 + x9
 
Fit Measure
─────────────────────────────────────────────────────────────────────────────────────
      Χ²      DF          P    CFI  RMSEA   SRMR    TLI       AIC       BIC      BIC2
─────────────────────────────────────────────────────────────────────────────────────
  85.306  24.000  0.000 ***  0.931  0.092  0.065  0.896  7517.490  7595.339  7528.739
─────────────────────────────────────────────────────────────────────────────────────
*** p < 0.001, ** p < 0.01, * p < 0.05, + p < 0.1

 
Factor Loadings
────────────────────────────────────────────────────────────────────────────────
  Latent.Factor  Observed.Var  Std.Est     SE       Z          P          95% CI
────────────────────────────────────────────────────────────────────────────────
            DV1            x1    0.772  0.055  14.041  0.000 ***  [0.664, 0.880]
                           x2    0.424  0.060   7.105  0.000 ***  [0.307, 0.540]
                           x3    0.581  0.055  10.539  0.000 ***  [0.473, 0.689]
            DV2            x4    0.852  0.023  37.776  0.000 ***  [0.807, 0.896]
                           x5    0.855  0.022  38.273  0.000 ***  [0.811, 0.899]
                           x6    0.838  0.023  35.881  0.000 ***  [0.792, 0.884]
            DV3            x7    0.570  0.053  10.714  0.000 ***  [0.465, 0.674]
                           x8    0.723  0.051  14.309  0.000 ***  [0.624, 0.822]
                           x9    0.665  0.051  13.015  0.000 ***  [0.565, 0.765]
────────────────────────────────────────────────────────────────────────────────
*** p < 0.001, ** p < 0.01, * p < 0.05, + p < 0.1

 
Model Covariances
──────────────────────────────────────────────────────────────
  Var.1  Var.2    Est     SE      Z          P          95% CI
──────────────────────────────────────────────────────────────
    DV1    DV2  0.459  0.064  7.189  0.000 ***  [0.334, 0.584]
    DV1    DV3  0.471  0.073  6.461  0.000 ***  [0.328, 0.613]
    DV2    DV3  0.283  0.069  4.117  0.000 ***  [0.148, 0.418]
──────────────────────────────────────────────────────────────
*** p < 0.001, ** p < 0.01, * p < 0.05, + p < 0.1

 
Model Variance
──────────────────────────────────────────────────────
  Var    Est     SE       Z          P          95% CI
──────────────────────────────────────────────────────
   x1  0.404  0.085   4.763  0.000 ***  [0.238, 0.571]
   x2  0.821  0.051  16.246  0.000 ***  [0.722, 0.920]
   x3  0.662  0.064  10.334  0.000 ***  [0.537, 0.788]
   x4  0.275  0.038   7.157  0.000 ***  [0.200, 0.350]
   x5  0.269  0.038   7.037  0.000 ***  [0.194, 0.344]
   x6  0.298  0.039   7.606  0.000 ***  [0.221, 0.374]
   x7  0.676  0.061  11.160  0.000 ***  [0.557, 0.794]
   x8  0.477  0.073   6.531  0.000 ***  [0.334, 0.620]
   x9  0.558  0.068   8.208  0.000 ***  [0.425, 0.691]
  DV1  1.000  0.000     NaN    NaN      [1.000, 1.000]
  DV2  1.000  0.000     NaN    NaN      [1.000, 1.000]
  DV3  1.000  0.000     NaN    NaN      [1.000, 1.000]
──────────────────────────────────────────────────────
*** p < 0.001, ** p < 0.01, * p < 0.05, + p < 0.1

 
Goodness of Fit:
 Warning. Poor χ² fit (p < 0.05). It is common to get p < 0.05. Check other fit measure.
 OK. Acceptable CFI fit (CFI > 0.90)
 Warning. Poor RMSEA fit (RMSEA > 0.08)
 OK. Good SRMR fit (SRMR < 0.08)
 Warning. Poor TLI fit (TLI < 0.90)
 OK. Barely acceptable factor loadings (0.4 < some loadings < 0.7)

Knit to R Markdown

if you want to produce these beautiful output in R Markdown. Calls this function and see the most up-to-date advice.

knit_to_Rmd()
OK. Required package "fansi" is installed

Note: To knit Rmd to HTML, add the following line to the setup chunk of your Rmd file: 
 "old.hooks <- fansi::set_knit_hooks(knitr::knit_hooks)"

Note: Use html_to_pdf to convert it to PDF. See ?html_to_pdf for more info

Ending

This conclude my briefed discussion of this package. I hope you enjoy the package, and please let me know if you have any feedback. If you like it, please considering giving a star on GitHub. Thank you.