xiaoruizhu/passo

An implementation of the unified framework for assessing Partial Association between ordinal variables after adjusting for a set of covariates (Dungang Liu, Shaobo Li, Yan Yu and Irini Moustaki (2020), accepted by the Journal of the American Statistical Association). This package provides a set of tools to quantify, visualize, and test partial associations between multiple ordinal variables.

PAsso: an R Package for Assessing Partial Association between Ordinal Variables

Overview

An implementation of the unified framework for assessing Parrtial Association between ordinal variables after adjusting for a set of covariates (Dungang Liu, Shaobo Li, Yan Yu and Irini Moustaki (2020), accepted by the Journal of the American Statistical Association). This package provides a set of tools to quantify, visualize, and test partial associations between multiple ordinal variables. It can produce a number of $\phi$ measures, partial regression plots, 3-D plots, and $p$-values for testing $H_0: \phi=0$ or $H_0: \phi \leq \delta$

Installation

The PAsso package is currently available on PAsso CRAN.

Install PAsso development version from GitHub (recommended)

# Install the development version from GitHub
if (!requireNamespace("devtools")) install.packages("devtools")
devtools::install_github("XiaoruiZhu/PAsso")

Install PAsso from the CRAN

# Install from CRAN
install.packages("PAsso")

# For macOS, if you have "error: 'math.h' file not found" for installing PAsso v0.1.9,
# the solution could be:
install.packages("https://cran.r-project.org/bin/macosx/el-capitan/contrib/3.6/PAsso_0.1.9.tgz", 
                  repos = NULL, type = "source")
                  
# If error "there is no package called 'gsl'" comes, try:
install.packages("gsl", type = "mac.binary")

Example

The following example shows the R code for evaluating the partial association between a binary variable $\text{PreVote.num}$ and a ordinal variable $\text{PID}$, while adjusting for age, education, and income. Specifically, $\text{PreVote.num}$ is the respondent's voting preference between Donald Trump and Hilary Clinton. And $\text{PID}$ is the respondent’s party identification with 7 ordinal levels from strong democrat (=1) to strong republication (=7). The data set is drawn from the 2016 American National Election Study.

library(PAsso)

PAsso_1 <- PAsso(responses = c("PreVote.num", "PID"),
                 adjustments = c("income.num", "age", "edu.year"),
                 data = ANES2016,
                 uni.model = "probit",
                 method = c("kendall"))

# Print the partial association matrix only
print(PAsso_1, 5)

# Provide:
# 1. partial association matrix;
# 2. marginal association matrix for comparison purpose;
# 3. summary of models' coefficients for model diagnostics and interpretation
summary(PAsso_1, 4)

# Plot partial association regression plot: residuals
plot(PAsso_1)

# Retrieve residuals that are used as ingredients for partial assocaition analyses
test_resids <- residuals(PAsso_1, draw = 1)
head(test_resids)

# test function: Conduct inference based on object of "PAsso.test" class ----------------------------
library(progress);

system.time(Pcor_SR_test1 <- test(object = PAsso_1, bootstrap_rep = 100, H0 = 0, parallel = F))
print(Pcor_SR_test1, digits=6)

# diagnostic.plot function -----------------------------------------------------
check_qq <- diagnostic.plot(object = PAsso_1, output = "qq")

check_fitted <- diagnostic.plot(object = PAsso_1, output = "fitted")

check_covar <- diagnostic.plot(object = PAsso_1, output = "covariate")

# Or more specific, draw residual-vs-covariate plot for the second model with
# response "PID" and covariate "income.num" 
diagnostic.plot(object = PAsso_1, output = "covariate", x_name = "income.num", model_id = 2)

# general association measure and 3-D plot for VOTE and PID ------------------
library("copula"); library("plotly")

# Draw all pairs
testPlots <- plot3D(PAsso_1)
testPlots$`PreVote.num v.s. PID`

# Draw just one pair
testPlots2 <- plot3D(object = PAsso_1, y1 = "PreVote.num", y2 = "PID")
testPlots2

# "PAsso" advanced using of the function: Input a few models directly ------------------------------
fit_vote <- glm(PreVote.num ~ income.num + age + edu.year, data = ANES2016,
                family = binomial(link = "probit"))
summary(fit_vote)

library(MASS)
fit_PID <- polr(as.factor(PID) ~ income.num + age + edu.year, data = ANES2016,
                method = "probit", Hess = TRUE)

summary(fit_PID)

system.time(PAsso_adv1 <- PAsso(fitted.models=list(fit_vote, fit_PID),
                                association = c("partial"),
                                method = c("kendall"),
                                resids.type = "surrogate")
)

# Partial association coefficients 
print(PAsso_adv1, digits = 3)
summary(PAsso_adv1, digits = 3)

References

Li et al. (2021). PAsso: an R Package for Assessing Partial Association between Ordinal Variables. The R Journal, 13(2), 239--252, <10.32614/RJ-2021-088>

Liu, D., Li, S., Yu, Y., & Moustaki, I. (2020). Assessing partial association between ordinal variables: quantification, visualization, and hypothesis testing. Journal of the American Statistical Association, 1-14. <10.1080/01621459.2020.1796394>

Liu, D., & Zhang, H. (2018). Residuals and diagnostics for ordinal regression models: A surrogate approach. Journal of the American Statistical Association, 113(522), 845-854. <10.1080/01621459.2017.1292915>

Greenwell, B.M., McCarthy, A.J., Boehmke, B.C. & Liu, D. (2018) Residuals and diagnostics for binary and ordinal regression models: An introduction to the sure package. The R Journal. <10.32614/RJ-2018-004>

passo's Issues

Avoid depending on the entire tidyverse package

Got suggestion from Hadley about this. I need to carefully read the Advanced R again. LOL

corr function for marginal correlation shouldn't need adjustments

When running the example in README, if use corr function for marginal correlation, it throws into error due to missing of "adjustments".

Fix it to not require the adjustments if run "marginal" association.

PAsso::residuals() for "glm" class object with "binomial" family will not generate surrogate residual by default

By default, we want the "residuals()" function in the PAsso package can include the features of the function "stats::residuals()". But it seems the default type of residuals generated from "PAsso::residuals()" for the "glm" class object with the "binomial" family is "deviance" residuals, rather than the surrogate residuals. This needs to be fixed.

residuals()

data("ANES2016")
fit.PID1<- polr(as.factor(PID)~age+edu.year+income.num, data=ANES2016, method="probit")
fit.PID2<- vglm(as.numeric(PID)~age+edu.year+income.num, data=ANES2016, family=acat(reverse=TRUE, parallel=TRUE))
r1<- residuals(fit.PID2, jitter="uniform",jitter.uniform.scale="response", nsim=30)
r2<- residuals(fit.PID2, jitter="uniform",jitter.uniform.scale="probability", nsim=30)
attr(r1,"draws")[,1]

the use of attr() is extremely inconvenient to get a simple output. Suggest to use more traditional way, e.g., r1[,1]
5th line, r2<-... does not work even if fit.PID2 is adjacent category logit model...

residuals()

Am I too annoying? hahaha~

a1<- residuals(fit.PID1, jitter="latent",jitter.uniform.scale="response")
a2<- residuals(fit.PID1, jitter="uniform",jitter.uniform.scale="response")
a3<- residuals(fit.PID1, jitter="latent",jitter.uniform.scale="probability")
a4<- residuals(fit.PID1, jitter="uniform",jitter.uniform.scale="probability")
hist(a1)
hist(a2)
hist(a3)
hist(a4)

a4 does not work. a1 and a3 seem to be the same thing, but what is a2???

Release PAsso 0.1.11

Prepare for release:

• git pull
• Check current CRAN check results
• Polish NEWS
• urlchecker::url_check()
• devtools::check(remote = TRUE, manual = TRUE)
• devtools::check_win_devel()
• rhub::check_for_cran()
• rhub::check(platform = 'ubuntu-rchk')
• rhub::check_with_sanitizers()
• revdepcheck::revdep_check(num_workers = 4)
• Update cran-comments.md
• git push

Submit to CRAN:

• usethis::use_version('patch')
• devtools::submit_cran()
• Approve email

Wait for CRAN...

• Accepted 🎉
• git push
• usethis::use_github_release()
• usethis::use_dev_version()
• git push

Install issue on macOS Mojave with R 3.6.3

As a maintainer, I test PAsso v0.1.8 on my macOS Mojave 10.14.5 with R 3.6.3. It throws issue when install:

fatal error: 'math.h' file not found

Then, the solution will be:

install.packages("https://cran.r-project.org/bin/macosx/el-capitan/contrib/3.6/PAsso_0.1.8.tgz", repos = NULL, type = "source")

However, on my end, another error comes:

there is no package called 'gsl'

Then, use binary installation of gsl package will solve the problem.

install.packages("gsl", type = "mac.binary")

Issue when output = "covariate" for diagnostic.plot()

data("ANES2016")
# Partial association between "PreVote.num" and "PID" after adjusting "income.num", "age", "edu.year"
# "PAsso" function: Only need input responses, adjustments, data
PAsso_1 <- PAsso(responses = c("PreVote.num", "PID"),
adjustments = c("income.num", "age", "edu.year"),
data = ANES2016,
uni.model = "probit",
method = c("kendall"))
diagnostic.plot(object = PAsso_1, output = "covariate", x_name = "income.num", model_id = 2)

Error in assign(x = x_name, object$data[, x_name]) : invalid first argument

diagnostic.plot() for adjacent category logit model

phi2 <- PAsso(responses = c("PID", "selfLR", "TrumpLR", "ClinLR"), adjustments = c("age", "edu.year", "income.num"), data = ANES2016, method = "kendall", model = c("acat", "acat", "acat", "acat"))
diagnostic.plot(phi2, output = "fitted")

The residuals are supposed to be in U(0,1) scale under acat model, but the plot shows different scale.