pnormmix <- function(q, mean1 = 0, sd1 = 1, mean2 = 0, sd2 = 1, pmix1 = 0.5,
lower.tail = TRUE) {
# Distribution function (pdf) of a mixture of two normal distributions.
#
# Args:
# q: vector of quantiles.
# mean1: mean of first normal distribution.
# sd1: standard deviation of first normal distribution.
# mean2: mean of second normal distribution.
# sd2: standard deviation of second normal distribution.
# pmix1: mixing proportion for the first distribution. Should be a
# number in the range (0, 1).
# lower.tail: logical; if TRUE(default), probabilities are P[X <= x];
# otherwise, P[X > x].
#
# Returns:
# Cumulative probability of q or 1 minus it on the mixture normal
# distribution.
# Error handling
stopifnot(pmix1 > 0, pmix1 < 1)
as.vector(c(pmix1, 1 - pmix1) %*%
sapply(q, pnorm, mean = c(mean1, mean2), sd = c(sd1, sd2),
lower.tail = lower.tail))
}
qnormmix <- function(p, mean1 = 0, sd1 = 1, mean2 = 0, sd2 = 1, pmix1 = 0.5,
lower.tail = TRUE) {
# Quantile function of a mixture of two normal distributions.
#
# Args:
# p: vector of probabilities.
# mean1: mean of first normal distribution.
# sd1: standard deviation of first normal distribution.
# mean2: mean of second normal distribution.
# sd2: standard deviation of second normal distribution.
# pmix1: mixing proportion for the first distribution. Should be a
# number in the range (0, 1).
# lower.tail: logical; if TRUE(default), probabilities are P[X <= x];
# otherwise, P[X > x].
#
# Returns:
# Quantile corresponding to p or 1 - p on the mixture normal
# distribution.
# Error handling
stopifnot(pmix1 > 0, pmix1 < 1, p >= 0, p <= 1)
f <- function(x) (pnormmix(x, mean1, sd1, mean2, sd2, pmix1,
lower.tail) - p)^2
start <- as.vector(c(pmix1, 1 - pmix1) %*%
sapply(p, qnorm, c(mean1, mean2), c(sd1, sd2),
lower.tail = lower.tail))
nlminb(start, f)$par
}
.bvnorm_kernel <- function(x, y, mu_x = 0, mu_y = 0, sd_x = 1, sd_y = 1,
cov_xy = 0) {
# Helper funcction for computing the kernel for bivariate normal density
cor <- cov_xy / sd_x / sd_y
numer <- (x - mu_x)^2 / sd_x^2 + (y - mu_y)^2 / sd_y^2 -
2 * cor * (x - mu_x) * (y - mu_y) / sd_x / sd_y
numer / (1 - cor^2)
}
contour_bvnorm <- function(mean1 = 0, sd1 = 1, mean2 = 0, sd2 = 1,
cor12 = 0, cov12 = NULL,
density = .95, length_out = 101,
bty = "L",
...) {
# Plot contour for a bivariate normal distribution
#
# Args:
# mean1: mean of first normal distribution (on x-axis).
# sd1: standard deviation of first normal distribution.
# mean2: mean of second normal distribution (on y-axis).
# sd2: standard deviation of second normal distribution.
# cor12: correlation in the bivariate normal.
# cov12: covariance in the bivariate normal. If not input, compute the
# covariance using the correlation and the standard deviations.
# density: density level, i.e., probability enclosed by the ellipse.
# length_out: number of values on the x-axis and on the y-axis to be
# evaluated; default to 101.
# bty: argument passed to the `contour` function.
# ...: other arguments passed to the `countour` funcction.
#
# Returns:
# a plot showing the contour of the bivariate normal distribution on
# a two-dimensional space.
# Error handling
stopifnot(cor12 >= -1, cor12 <= 1)
if (is.null(cov12)) cov12 <- cor12 * sd1 * sd2s
x_seq <- mean1 + seq(-3, 3, length.out = length_out) * sd1
y_seq <- mean2 + seq(-3, 3, length.out = length_out) * sd2
z <- outer(x_seq, y_seq, .bvnorm_kernel, mu_x = mean1, mu_y = mean2,
sd_x = sd1, sd_y = sd2, cov_xy = cov12)
contour(x_seq, y_seq, z, levels = qchisq(density, 2), drawlabels = FALSE,
bty = bty, ...)
}
.partit_bvnorm <- function(cut1, cut2, mean1 = 0, sd1 = 1, mean2 = 0, sd2 = 1,
cor12 = 0, cov12 = cor12 * sd1 * sd2) {
# Helper funcction for computing summary statistics from a selection approach
Sigma <- matrix(c(sd1^2, cov12, cov12, sd2^2), nrow = 2)
C <- pmnorm(c(cut1, cut2), c(mean1, mean2), Sigma)
B <- pnorm(cut1, mean1, sd1) - C
D <- pnorm(cut2, mean2, sd2) - C
A <- 1 - B - C - D
propsel <- A + B
success_ratio <- A / propsel
sensitivity <- A / (A + D)
specificity <- C / (C + B)
c(A, B, C, D, propsel, success_ratio, sensitivity, specificity)
}
PartInv <- function(propsel, kappa_r, kappa_f = kappa_r,
phi_r, phi_f = phi_r, lambda_r, lambda_f = lambda_r,
Theta_r, Theta_f = Theta_r, tau_r, tau_f = tau_r,
pmix_ref = 0.5, plot_contour = TRUE, ...) {
# Evaluate partial measurement invariance using Millsap & Kwok's (2004)
# approach
#
# Args:
# propsel: proportion of selection.
# kappa_r: latent factor mean for the reference group.
# kappa_f: (optional) latent factor mean for the focal group;
# if no input, set equal to kappa_r.
# phi_r: latent factor variance for the reference group.
# phi_f: (optional) latent factor variance for the focal group;
# if no input, set equal to phi_r.
# lambda_r: a vector of factor loadings for the reference group.
# lambda_f: (optional) a vector of factor loadings for the focal group;
# if no input, set equal to lambda_r.
# tau_r: a vector of measurement intercepts for the reference group.
# tau_f: (optional) a vector of measurement intercepts for the focal group;
# if no input, set equal to tau_r.
# Theta_r: a matrix of the unique factor variances and covariances
# for the reference group.
# Theta_f: (optional) a matrix of the unique factor variances and
# covariances for the focal group;
# if no input, set equal to Theta_r.
# pmix_ref: Proportion of the reference group;
# default to 0.5 (i.e., two populations have equal size).
# plot_contour: logical; whether the contour of the two populations
# should be plotted; default to TRUE.
# ...: other arguments passed to the `countour` funcction.
#
# Returns:
# a list of four elements and a plot if plot_contour == TRUE:
# - propsel: echo the same argument as input.
# - cutpt_xi: cut point on the latent scale (xi).
# - cutpt_z: cut point on the observed scale (Z).
# - summary: A 8 x 2 table, with columns representing the reference
# and the focal groups, and the rows represent
# probabilities of true positive (A), false positive (B),
# true negative (C), false negative (D); proportion selected,
# success ratio, sensitivity, and specificity.
# Error handling
stopifnot(length(kappa_r) == 1, length(kappa_f) == 1, length(phi_r) == 1,
length(phi_f) == 1)
if (length(Theta_r) == length(lambda_r)) Theta_r <- diag(Theta_r)
if (length(Theta_f) == length(lambda_f)) Theta_f <- diag(Theta_f)
library(mnormt) # load `mnormt` package
mean_zr <- sum(tau_r) + sum(lambda_r) * kappa_r
mean_zf <- sum(tau_f) + sum(lambda_f) * kappa_f
sd_zr <- sqrt(sum(lambda_r)^2 * phi_r + sum(Theta_r))
sd_zf <- sqrt(sum(lambda_f)^2 * phi_f + sum(Theta_f))
cov_z_xir <- sum(lambda_r) * phi_r
cov_z_xif <- sum(lambda_f) * phi_f
sd_xir <- sqrt(phi_r)
sd_xif <- sqrt(phi_f)
cut_xi <- qnormmix(propsel, kappa_r, sd_xir, kappa_f, sd_xif,
pmix_ref, lower.tail = FALSE)
cut_z <- qnormmix(propsel, mean_zr, sd_zr, mean_zf, sd_zf,
pmix_ref, lower.tail = FALSE)
partit_1 <- .partit_bvnorm(cut_xi, cut_z, kappa_r, sd_xir, mean_zr, sd_zr,
cov12 = cov_z_xir)
partit_2 <- .partit_bvnorm(cut_xi, cut_z, kappa_f, sd_xif, mean_zf, sd_zf,
cov12 = cov_z_xif)
dat <- data.frame("Reference" = partit_1, "Focal" = partit_2,
row.names = c("A (true positive)", "B (false positive)",
"C (true negative)", "D (false negative)",
"Proportion selected", "Success ratio",
"Sensitivity", "Specificity"))
if (plot_contour) {
x_lim <- range(c(kappa_r + c(-3, 3) * sd_xir,
kappa_f + c(-3, 3) * sd_xif))
y_lim <- range(c(mean_zr + c(-3, 3) * sd_zr,
mean_zf + c(-3, 3) * sd_zf))
contour_bvnorm(kappa_r, sd_xir, mean_zr, sd_zr, cov12 = cov_z_xir,
xlab = bquote("Latent Score" ~ (xi)),
ylab = bquote("Observed Composite" ~ (italic(Z))),
lwd = 2, col = "red", xlim = x_lim, ylim = y_lim,
...)
contour_bvnorm(kappa_f, sd_xif, mean_zf, sd_zf, cov12 = cov_z_xif,
add = TRUE, lty = "dashed", lwd = 2, col = "blue",
...)
legend("topleft", c("Reference group", "Focal group"),
lty = c("solid", "dashed"), col = c("red", "blue"))
abline(h = cut_z, v = cut_xi)
x_cord <- rep(cut_xi + c(.25, -.25) * sd_xir, 2)
y_cord <- rep(cut_z + c(.25, -.25) * sd_zr, each = 2)
text(x_cord, y_cord, c("A", "B", "D", "C"))
}
list(propsel = propsel, cutpt_xi = cut_xi, cutpt_z = cut_z,
summary = round(dat, 3))
}
##########
mytest<-function(propsel, kappa_r, kappa_f, phi_r, phi_f,
lambda_r,
tau_r,
tau_f,
Theta_r,
Theta_f,
pmix_ref){PartInv(propsel, kappa_r, kappa_f, phi_r, phi_f,
lambda_r,
tau_r,
tau_f,
Theta_r,
Theta_f,
pmix_ref)}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.