混合Cor（mydat，method =“pearson”）中的错误：在R中

Question

我尝试执行点双序列相关。我使用mixedCor 来自library("psych")

library("psych")
mixedCor(mydata,method="pearson")

我得到了错误

mixedCor(mydat, method = "pearson") 中的错误： 我试图弄清楚哪些是连续的，哪些是多分的，但失败了。请指定 x、p 和 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= ---

Answer 1

您需要在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