---
title: "Introduction to the pervasive package"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to the pervasive package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```


The pervasive package is designed to help assess the pervasiveness of effects in correlational data by providing the Observed Percentage of Concordant Pairs (OPCP) and association rule mining for binned variables. 

If researchers want OPCP values without conducting association rule mining, they may use  OPCP(), OPCP_glm() for binary outcomes, or OPCP_mat() for a whole matrix. These OPCP functions do not apply any binning.

The pervasive_tric(), pervasive_tric_glm(), pervasive_dic() and pervasive_dic_glm() all provide the OPCP and the results of targeted association rule mining. The function pervasive_tric() should likely be considered the go-to for most research contexts. When the sample size is too small for trichotomization to make sense, or if too many variables are binary, researchers may use pervasive_dic() instead. The functions pervasive_tric_glm() and pervasive_dic_glm() should be preferred if the outcome is binary. 

This vignette will walk you through examples for the pervasive functions using data from the psychTools package.


```{r}
#> Install from CRAN (not yet available)
#install.packages("pervasive")

#> Install the development version from GitHub
#devtools::install_github("dlajoiemoncton/pervasive")

#load the package
library(pervasive)

```

We can start with a quick overview of the OPCP_mat(), OPCP(), and OPCP_glm() functions. OPCP is pretty easy to calculate for any models (correlate predicted values and observed values using the "kendall" method, divide the result by 2 and add .5), but for convenience, the pervasive package provides functions for linear and logistic regression. In this vignette, we will use data from the psychTools package to illustrate the functions.  


```{r, message=FALSE}
#>Example using the spi dataset from the psychTools package. 

#>Big 5 factor scores are calculated as suggested in the psych package documentation. Variables should be scored prior to the analysis.
sc <- psych::scoreVeryFast(psychTools::spi.keys, psychTools::spi)

#>scores for traits are combined to the original dataset
spi_sc <- cbind(psychTools::spi, sc)

#>Let's use age as the outcome of interest and the Big 5 as predictors for the regression models and sex as the outcome for binomial logistic regression models

spi_sc_age_sex_B5 <- spi_sc |>
 dplyr::select(age, sex, Agree, Consc, Neuro, Extra, Open) |> 
  na.omit()

spi_sc_age <- spi_sc_age_sex_B5 |>
dplyr::select(age, Agree, Consc, Neuro, Extra, Open)

spi_sc_sex <- spi_sc_age_sex_B5 |>
 dplyr::select(sex, Agree, Consc, Neuro, Extra, Open)
 
#>The glm() function expects outputs to be coded as 0 and 1 but sex is coded 1 and 2. There are some NAs for sex.
spi_sc_sex$sex = spi_sc_sex$sex -1
```



The OPCP_mat function works like cor(). Instead of returning a correlation matrix, it provides a matrix containing Pearson's r in the lower triangle and OPCP values in the upper triangle. It takes a little longer to run than a normal correlation table. This may be worth considering if you have a quite large number of variables. In the resulting matrix we can see, for example, that in approcximately 56% of pairs of participants (excluding ties), the person with the higher agreeableness score is also the older participant. 

```{r}
OPCP_mat(spi_sc_age_sex_B5)


```


The other pervasive functions expect a formula parameter that is in the same shape as a typical lm() or glm() model. 

```{r}
formula <- age ~ Agree + Consc + Neuro + Extra + Open 

formula_glm <- sex ~ Agree + Consc + Neuro + Extra + Open

#>This would also be acceptable: formula <- formula(age ~ Agree + Consc + Neuro + Extra + Open)
#>It is possible to include a single predictor if one is working in a bivariate case 

```


OPCP() and OPCP_glm() can provide OPCPs for regression models. In the context of binomial logistic regression, researchers have options for model performance evaluation based on the confusion matrix (e.g., accuracy, specificity, Jaccard index, etc.) such that the OPCP might feel like simply one of many. One advantage for the OPCP is that it makes it convenient to directly compare the performance of glm() and lm() models, on the same metric. In the examples below, there appears to be more information in personality for sex than for age.

```{r}
OPCP(formula = formula, data = spi_sc_age)
OPCP_glm(formula = formula_glm, data = spi_sc_sex)
```

The following functions (pervasive_tric(), _dic(), etc.) include some results from association rule mining. This allows some data exploration and provides context for models. A new parameter is required by the functions, min_support.

This parameter refers to the minimum frequency with which a configuration appears in the data. The idea is to go low enough where we do not miss interesting patterns but high enough that configurations are not completely spurious. 

If a dataset is smaller (everything is contextual, but, say, below 300), it may be worthwhile to consider using pervasive_dic() instead of pervasive_tric() and to use a higer min_support. The old rule-of-thumb of n > 30 for the smallest class is an easy "minimum" to aim for, but context always matters. Here, .03*3946 = 118.38, such that configurations with fewer than 119 individuals will not be considered to have enough support to be of note. 

```{r, message=FALSE}
example <- pervasive_tric(formula = formula, data = spi_sc_age, min_support = .03)

example


#> or print(pervasive_tric(formula = formula, data = spi_sc_age, min_support = .03))
```

From the results, it appears we would be rather unlikely to meet individuals with the patterns of personality traits suggested for old and young people by a linear regression when data is trichotomized. The regression model explains ~9% of variance. This coincides with an observed percentage of concordant pairs of ~61%, as reported above. The middle category of age appears to be rather difficult to predict, with a lift of 1.20 being the best that is achieved (lift can globally be interpreted in the same way as an overall or omnibus odds ratio). 

In the frequency tables below, we can see that the bins are of fairly equivalent size. However, the cutoffs make for some unequal ranges for many variables. This is typical with skewed variables.  

```{r}
example$freq_tables
```


We can compare to the results of dichotomizing instead of trichotomizing. This improves support for the regression solutions. However, lift is generally smaller.

```{r, message=FALSE}

example_dic <- pervasive_dic(formula = formula, data = spi_sc_age, min_support = .03)

example_dic

example_dic$freq_tables

```



If we want to examine pervasiveness with binary outcomes, we can use pervasive_tric_glm() or pervasive_dic_glm().

When variables are trichotomized, we can see that the "prototypical" woman suggested by the regression, with high agreeableness, conscienciousness, neuroticism, extraversion and low openness, is again rather rare (14 individuals out of a sample of 3946). This does not imply that the linear model is somehow "wrong", it simply provides context to interpret the results.  

```{r, message=FALSE}

 example_dic_glm <- pervasive_dic_glm(formula = formula_glm, data = spi_sc_sex, min_support = .03)
 example_dic_glm
 
 example_tric_glm <- pervasive_tric_glm(formula = formula_glm, data = spi_sc_sex, min_support = .03)
 example_tric_glm
```



