mixedCor(mydat, method = "pearson") 错误：在 R 中
Error in mixedCor(mydat, method = "pearson") : in R
我尝试执行点双序列相关。我使用 mixedCor 来自 library("psych")
library("psych")
mixedCor(mydata,method="pearson")
我得到了错误
Error in mixedCor(mydat, method = "pearson") :
I tried to figure out which where continuous and which were polytomous, but failed. Please try again by specifying x, p, and d.
数据
mydata(dput)=structure(list(x1 = c(9L, 7L, 2L, 5L, 6L, 5L, 8L, 2L, 4L, 5L,
8L, 3L, 2L, 3L, 2L, 5L, 4L, 5L, 4L, 4L, 4L), x2 = c(6L, 1L, 3L,
1L, 7L, 5L, 3L, 3L, 2L, 2L, 3L, 4L, 5L, 3L, 9L, 5L, 6L, 4L, 3L,
6L, 7L), x3 = c(8L, 7L, 7L, 9L, 8L, 6L, 7L, 7L, 9L, 7L, 9L, 7L,
11L, 9L, 7L, 10L, 3L, 10L, 8L, 6L, 6L), y = c(0L, 0L, 1L, 0L,
0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L,
1L)), .Names = c("x1", "x2", "x3", "y"), class = "data.frame", row.names = c(NA,
-21L))
x1,x2,x3 是比例变量，y 是分类变量（0 或 1）
如何执行这种类型的相关，输出将具有这样的 p 值？
y x3 x2 x1
y 1,0000 ,1114 ,2201 -,2597
p= --- p=,631 p=,338 p=,256
x3 ,1114 1,0000 -,2630 ,0457
p=,631 p= --- p=,249 p=,844
x2 ,2201 -,2630 1,0000 -,1548
p=,338 p=,249 p= --- p=,503
x1 -,2597 ,0457 -,1548 1,0000
p=,256 p=,844 p=,503 p= ---
您需要在 mixedCor:
中指定 c 和 d（分别为连续变量集和二分变量集）
library(psych)
mixedCor(data=mydata, c=1:3, d=4, method="pearson")
# x1 x2 x3 y
# x1 1.00
# x2 -0.13 1.00
# x3 0.03 -0.25 1.00
# y -0.32 0.26 0.07 1.00
如果需要计算p值，可以使用下面的mycor.ci函数（cor.ci的修改版）：
mycor.ci <- function (x, keys = NULL, n.iter = 100, p = 0.05, overlap = FALSE,
poly = FALSE, method = "pearson", plot = TRUE, minlength = 5,
cvars=NULL, pvars=NULL, dvars=NULL, ...) {
cl <- match.call()
n.obs <- dim(x)[1]
if (is.null(keys) && overlap)
overlap <- FALSE
if (poly) {
ncat <- 8
nvar <- dim(x)[2]
tx <- table(as.matrix(x))
if (dim(tx)[1] == 2) {
tet <- tetrachoric(x)
typ = "tet"
}
else {
tab <- apply(x, 2, function(x) table(x))
if (is.list(tab)) {
len <- lapply(tab, function(x) length(x))
}
else {
len <- dim(tab)[1]
}
if (is.null(dvars)) dvars <- subset(1:nvar, len == 2)
if (is.null(pvars)) pvars <- subset(1:nvar, ((len > 2) & (len <= ncat)))
if (is.null(cvars)) cvars <- subset(1:nvar, (len > ncat))
if (length(pvars) == ncol(x)) {
tet <- polychoric(x)
typ = "poly"
}
else {
plot(pvars)
tet <- mixedCor(data=x, c=cvars, d=dvars, method="pearson")
typ = "mixed"
}
}
Rho <- tet$rho
}
else {
Rho <- cor(x, use = "pairwise", method = method)
}
if (!is.null(keys)) {
bad <- FALSE
if (any(is.na(Rho))) {
warning(sum(is.na(Rho)), " of the item correlations are NA and thus finding scales that include those items will not work.\n We will try to do it for those scales which do not include those items.\n \n Proceed with caution. ")
bad <- TRUE
rho <- apply(keys, 2, function(x) colMeans(apply(keys,
2, function(x) colMeans(Rho * x, na.rm = TRUE)) *
x, na.rm = TRUE))
}
else {
rho <- t(keys) %*% Rho %*% keys
}
}
else {
rho <- Rho
}
if (overlap) {
key.var <- diag(t(keys) %*% keys)
var <- diag(rho)
n.keys <- ncol(keys)
key.av.r <- (var - key.var)/(key.var * (key.var - 1))
item.cov <- t(keys) %*% Rho
raw.cov <- item.cov %*% keys
adj.cov <- raw.cov
for (i in 1:(n.keys)) {
for (j in 1:i) {
adj.cov[i, j] <- adj.cov[j, i] <- raw.cov[i,
j] - sum(keys[, i] * keys[, j]) + sum(keys[,
i] * keys[, j] * sqrt(key.av.r[i] * key.av.r[j]))
}
}
diag(adj.cov) <- diag(raw.cov)
rho <- cov2cor(adj.cov)
}
rho <- cov2cor(rho)
nvar <- dim(rho)[2]
if (n.iter > 1) {
replicates <- list()
rep.rots <- list()
replicates <- parallel::mclapply(1:n.iter, function(XX) {
xs <- x[sample(n.obs, n.obs, replace = TRUE), ]
{
if (poly) {
if (typ != "tet") {
capture.output(tets <- mixedCor(data=xs, c=cvars, d=dvars, method="pearson"))
}
else {
tets <- tetrachoric(xs)
}
R <- tets$rho
}
else {
R <- cor(xs, use = "pairwise", method = method)
}
if (!is.null(keys)) {
if (bad) {
covariances <- apply(keys, 2, function(x) colMeans(apply(keys,
2, function(x) colMeans(R * x, na.rm = TRUE)) *
x, na.rm = TRUE))
}
else {
covariances <- t(keys) %*% R %*% keys
}
r <- cov2cor(covariances)
}
else {
r <- R
}
if (overlap) {
var <- diag(covariances)
item.cov <- t(keys) %*% R
raw.cov <- item.cov %*% keys
adj.cov <- raw.cov
key.av.r <- (var - key.var)/(key.var * (key.var -
1))
for (i in 1:(n.keys)) {
for (j in 1:i) {
adj.cov[i, j] <- adj.cov[j, i] <- raw.cov[i,
j] - sum(keys[, i] * keys[, j]) + sum(keys[,
i] * keys[, j] * sqrt(key.av.r[i] * key.av.r[j]))
}
}
diag(adj.cov) <- diag(raw.cov)
r <- cov2cor(adj.cov)
}
rep.rots <- r[lower.tri(r)]
}
})
}
replicates <- matrix(unlist(replicates), ncol = nvar * (nvar -
1)/2, byrow = TRUE)
z.rot <- fisherz(replicates)
means.rot <- colMeans(z.rot, na.rm = TRUE)
sds.rot <- apply(z.rot, 2, sd, na.rm = TRUE)
sds.rot <- fisherz2r(sds.rot)
ci.rot.lower <- means.rot + qnorm(p/2) * sds.rot
ci.rot.upper <- means.rot + qnorm(1 - p/2) * sds.rot
means.rot <- fisherz2r(means.rot)
ci.rot.lower <- fisherz2r(ci.rot.lower)
ci.rot.upper <- fisherz2r(ci.rot.upper)
low.e <- apply(replicates, 2, quantile, p/2, na.rm = TRUE)
up.e <- apply(replicates, 2, quantile, 1 - p/2, na.rm = TRUE)
tci <- abs(means.rot)/sds.rot
ptci <- pnorm(tci)
ci.rot <- data.frame(lower = ci.rot.lower, low.e = low.e,
upper = ci.rot.upper, up.e = up.e, p = 2 * (1 - ptci))
cnR <- abbreviate(colnames(rho), minlength = minlength)
k <- 1
for (i in 1:(nvar - 1)) {
for (j in (i + 1):nvar) {
rownames(ci.rot)[k] <- paste(cnR[i], cnR[j], sep = "-")
k <- k + 1
}
}
results <- list(rho = rho, means = means.rot, sds = sds.rot,
tci = tci, ptci = ptci, ci = ci.rot, Call = cl, replicates = replicates)
if (plot) {
cor.plot(rho, numbers = TRUE, cuts = c(0.001, 0.01, 0.05),
pval = 2 * (1 - ptci), ...)
}
class(results) <- c("psych", "cor.ci")
return(results)
}
语法是：
library(parallel)
set.seed(123)
mycor.ci(x=mydata, method="pearson", poly=TRUE, cvars=1:3, dvars=4, n.iter=1000)
# Coefficients and bootstrapped confidence intervals
# x1 x2 x3 y
# x1 1.00
# x2 -0.13 1.00
# x3 0.03 -0.25 1.00
# y -0.32 0.26 0.07 1.00
# scale correlations and bootstrapped confidence intervals
# lower.emp lower.norm estimate upper.norm upper.emp p
# x1-x2 -0.50 -0.51 -0.13 0.29 0.30 0.56
# x1-x3 -0.31 -0.31 0.03 0.40 0.38 0.79
# x1-y -0.79 -0.75 -0.32 0.24 0.19 0.26
# x2-x3 -0.58 -0.57 -0.25 0.12 0.13 0.19
# x2-y -0.34 -0.36 0.26 0.72 0.78 0.42
# x3-y -0.51 -0.48 0.07 0.57 0.56 0.84
我尝试执行点双序列相关。我使用 mixedCor 来自 library("psych")
library("psych")
mixedCor(mydata,method="pearson")
我得到了错误
Error in mixedCor(mydat, method = "pearson") : I tried to figure out which where continuous and which were polytomous, but failed. Please try again by specifying x, p, and d.
数据
mydata(dput)=structure(list(x1 = c(9L, 7L, 2L, 5L, 6L, 5L, 8L, 2L, 4L, 5L,
8L, 3L, 2L, 3L, 2L, 5L, 4L, 5L, 4L, 4L, 4L), x2 = c(6L, 1L, 3L,
1L, 7L, 5L, 3L, 3L, 2L, 2L, 3L, 4L, 5L, 3L, 9L, 5L, 6L, 4L, 3L,
6L, 7L), x3 = c(8L, 7L, 7L, 9L, 8L, 6L, 7L, 7L, 9L, 7L, 9L, 7L,
11L, 9L, 7L, 10L, 3L, 10L, 8L, 6L, 6L), y = c(0L, 0L, 1L, 0L,
0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L,
1L)), .Names = c("x1", "x2", "x3", "y"), class = "data.frame", row.names = c(NA,
-21L))
x1,x2,x3 是比例变量，y 是分类变量（0 或 1）
如何执行这种类型的相关，输出将具有这样的 p 值？
y x3 x2 x1
y 1,0000 ,1114 ,2201 -,2597
p= --- p=,631 p=,338 p=,256
x3 ,1114 1,0000 -,2630 ,0457
p=,631 p= --- p=,249 p=,844
x2 ,2201 -,2630 1,0000 -,1548
p=,338 p=,249 p= --- p=,503
x1 -,2597 ,0457 -,1548 1,0000
p=,256 p=,844 p=,503 p= ---
您需要在 mixedCor:
c 和 d（分别为连续变量集和二分变量集）
library(psych)
mixedCor(data=mydata, c=1:3, d=4, method="pearson")
# x1 x2 x3 y
# x1 1.00
# x2 -0.13 1.00
# x3 0.03 -0.25 1.00
# y -0.32 0.26 0.07 1.00
如果需要计算p值，可以使用下面的mycor.ci函数（cor.ci的修改版）：
mycor.ci <- function (x, keys = NULL, n.iter = 100, p = 0.05, overlap = FALSE,
poly = FALSE, method = "pearson", plot = TRUE, minlength = 5,
cvars=NULL, pvars=NULL, dvars=NULL, ...) {
cl <- match.call()
n.obs <- dim(x)[1]
if (is.null(keys) && overlap)
overlap <- FALSE
if (poly) {
ncat <- 8
nvar <- dim(x)[2]
tx <- table(as.matrix(x))
if (dim(tx)[1] == 2) {
tet <- tetrachoric(x)
typ = "tet"
}
else {
tab <- apply(x, 2, function(x) table(x))
if (is.list(tab)) {
len <- lapply(tab, function(x) length(x))
}
else {
len <- dim(tab)[1]
}
if (is.null(dvars)) dvars <- subset(1:nvar, len == 2)
if (is.null(pvars)) pvars <- subset(1:nvar, ((len > 2) & (len <= ncat)))
if (is.null(cvars)) cvars <- subset(1:nvar, (len > ncat))
if (length(pvars) == ncol(x)) {
tet <- polychoric(x)
typ = "poly"
}
else {
plot(pvars)
tet <- mixedCor(data=x, c=cvars, d=dvars, method="pearson")
typ = "mixed"
}
}
Rho <- tet$rho
}
else {
Rho <- cor(x, use = "pairwise", method = method)
}
if (!is.null(keys)) {
bad <- FALSE
if (any(is.na(Rho))) {
warning(sum(is.na(Rho)), " of the item correlations are NA and thus finding scales that include those items will not work.\n We will try to do it for those scales which do not include those items.\n \n Proceed with caution. ")
bad <- TRUE
rho <- apply(keys, 2, function(x) colMeans(apply(keys,
2, function(x) colMeans(Rho * x, na.rm = TRUE)) *
x, na.rm = TRUE))
}
else {
rho <- t(keys) %*% Rho %*% keys
}
}
else {
rho <- Rho
}
if (overlap) {
key.var <- diag(t(keys) %*% keys)
var <- diag(rho)
n.keys <- ncol(keys)
key.av.r <- (var - key.var)/(key.var * (key.var - 1))
item.cov <- t(keys) %*% Rho
raw.cov <- item.cov %*% keys
adj.cov <- raw.cov
for (i in 1:(n.keys)) {
for (j in 1:i) {
adj.cov[i, j] <- adj.cov[j, i] <- raw.cov[i,
j] - sum(keys[, i] * keys[, j]) + sum(keys[,
i] * keys[, j] * sqrt(key.av.r[i] * key.av.r[j]))
}
}
diag(adj.cov) <- diag(raw.cov)
rho <- cov2cor(adj.cov)
}
rho <- cov2cor(rho)
nvar <- dim(rho)[2]
if (n.iter > 1) {
replicates <- list()
rep.rots <- list()
replicates <- parallel::mclapply(1:n.iter, function(XX) {
xs <- x[sample(n.obs, n.obs, replace = TRUE), ]
{
if (poly) {
if (typ != "tet") {
capture.output(tets <- mixedCor(data=xs, c=cvars, d=dvars, method="pearson"))
}
else {
tets <- tetrachoric(xs)
}
R <- tets$rho
}
else {
R <- cor(xs, use = "pairwise", method = method)
}
if (!is.null(keys)) {
if (bad) {
covariances <- apply(keys, 2, function(x) colMeans(apply(keys,
2, function(x) colMeans(R * x, na.rm = TRUE)) *
x, na.rm = TRUE))
}
else {
covariances <- t(keys) %*% R %*% keys
}
r <- cov2cor(covariances)
}
else {
r <- R
}
if (overlap) {
var <- diag(covariances)
item.cov <- t(keys) %*% R
raw.cov <- item.cov %*% keys
adj.cov <- raw.cov
key.av.r <- (var - key.var)/(key.var * (key.var -
1))
for (i in 1:(n.keys)) {
for (j in 1:i) {
adj.cov[i, j] <- adj.cov[j, i] <- raw.cov[i,
j] - sum(keys[, i] * keys[, j]) + sum(keys[,
i] * keys[, j] * sqrt(key.av.r[i] * key.av.r[j]))
}
}
diag(adj.cov) <- diag(raw.cov)
r <- cov2cor(adj.cov)
}
rep.rots <- r[lower.tri(r)]
}
})
}
replicates <- matrix(unlist(replicates), ncol = nvar * (nvar -
1)/2, byrow = TRUE)
z.rot <- fisherz(replicates)
means.rot <- colMeans(z.rot, na.rm = TRUE)
sds.rot <- apply(z.rot, 2, sd, na.rm = TRUE)
sds.rot <- fisherz2r(sds.rot)
ci.rot.lower <- means.rot + qnorm(p/2) * sds.rot
ci.rot.upper <- means.rot + qnorm(1 - p/2) * sds.rot
means.rot <- fisherz2r(means.rot)
ci.rot.lower <- fisherz2r(ci.rot.lower)
ci.rot.upper <- fisherz2r(ci.rot.upper)
low.e <- apply(replicates, 2, quantile, p/2, na.rm = TRUE)
up.e <- apply(replicates, 2, quantile, 1 - p/2, na.rm = TRUE)
tci <- abs(means.rot)/sds.rot
ptci <- pnorm(tci)
ci.rot <- data.frame(lower = ci.rot.lower, low.e = low.e,
upper = ci.rot.upper, up.e = up.e, p = 2 * (1 - ptci))
cnR <- abbreviate(colnames(rho), minlength = minlength)
k <- 1
for (i in 1:(nvar - 1)) {
for (j in (i + 1):nvar) {
rownames(ci.rot)[k] <- paste(cnR[i], cnR[j], sep = "-")
k <- k + 1
}
}
results <- list(rho = rho, means = means.rot, sds = sds.rot,
tci = tci, ptci = ptci, ci = ci.rot, Call = cl, replicates = replicates)
if (plot) {
cor.plot(rho, numbers = TRUE, cuts = c(0.001, 0.01, 0.05),
pval = 2 * (1 - ptci), ...)
}
class(results) <- c("psych", "cor.ci")
return(results)
}
语法是：
library(parallel)
set.seed(123)
mycor.ci(x=mydata, method="pearson", poly=TRUE, cvars=1:3, dvars=4, n.iter=1000)
# Coefficients and bootstrapped confidence intervals
# x1 x2 x3 y
# x1 1.00
# x2 -0.13 1.00
# x3 0.03 -0.25 1.00
# y -0.32 0.26 0.07 1.00
# scale correlations and bootstrapped confidence intervals
# lower.emp lower.norm estimate upper.norm upper.emp p
# x1-x2 -0.50 -0.51 -0.13 0.29 0.30 0.56
# x1-x3 -0.31 -0.31 0.03 0.40 0.38 0.79
# x1-y -0.79 -0.75 -0.32 0.24 0.19 0.26
# x2-x3 -0.58 -0.57 -0.25 0.12 0.13 0.19
# x2-y -0.34 -0.36 0.26 0.72 0.78 0.42
# x3-y -0.51 -0.48 0.07 0.57 0.56 0.84