r - 多序列相关的 P 值

Question

polyserial()我有一些关于{polycor} 函数的基本问题。

1. rho 是否存在 p 值，还是可以计算出来？
2. 对于双变量正态的假设，测试的零假设是“是的，双变量正态”吗？也就是说，我想要一个高还是低的 p 值。

谢谢。

Answer 1

如果您使用以下方式形成返回的对象：

 polS <- polyserial(x, y, ML=TRUE, std.err=TRUE) # ML estimate

...您应该没有困难形成假设的 p 值：rho == 0使用由参数除以其标准误差的比率形成的 z 统计量。但这与检验双变量正态性假设不同。为此，您需要检查polS. 'polycor' 类对象的 print 方法用一个漂亮的小句子把它交给你。您以通常的方式解释该结果：低 p 值是反对零假设的有力证据（在本例中为 H0：双变量正态性）。作为一名科学家，你不“想要”任何一个结果。你想了解数据告诉你什么。

Answer 2

我给包作者发了电子邮件——因为我有同样的问题），根据他的澄清，我提供了我的答案：

首先，一个简单的问题：较高的 p 值（传统上 > 0.05）让您更有信心相信分布是双变量正态分布。较低的 p 值表示非正态分布，但是，如果样本量足够大，最大似然估计（选项ML=TRUE），非正态性无关紧要；无论如何，相关性仍然可靠。

现在，对于更难的问题：要计算 p 值，您需要：

1. 使用 std.err=TRUE 选项执行 polyserial 以访问更多详细信息。
2. 从生成的 polyserial 对象中，访问var[1, 1]元素。var是参数估计的协方差矩阵，并且sqrt(var[1, 1])是标准误差（在 rho 结果之后的输出中显示在括号中）。
3. 根据标准误差，您可以根据下面的 R 代码计算 p 值。

以下是一些代码，可根据多序列文档中的示例代码，用可复制的 R 代码说明这一点：

library(mvtnorm)
library(polycor)

set.seed(12345)
data <- rmvnorm(1000, c(0, 0), matrix(c(1, .5, .5, 1), 2, 2))
x <- data[,1]
y <- data[,2]
y <- cut(y, c(-Inf, -1, .5, 1.5, Inf))

# 2-step estimate
poly_2step <- polyserial(x, y, std.err=TRUE)  
poly_2step
## 
## Polyserial Correlation, 2-step est. = 0.5085 (0.02413)
## Test of bivariate normality: Chisquare = 8.604, df = 11, p = 0.6584
std.err_2step <- sqrt(poly_2step$var[1, 1])
std.err_2step
## [1] 0.02413489
p_value_2step <- 2 * pnorm(-abs(poly_2step$rho / std.err_2step))
p_value_2step
## [1] 1.529176e-98
# ML estimate
poly_ML <- polyserial(x, y, ML=TRUE, std.err=TRUE) 
poly_ML
## 
## Polyserial Correlation, ML est. = 0.5083 (0.02466)
## Test of bivariate normality: Chisquare = 8.548, df = 11, p = 0.6635
## 
##                  1      2       3
## Threshold -0.98560 0.4812 1.50700
## Std.Err.   0.04408 0.0379 0.05847
std.err_ML <- sqrt(poly_ML$var[1, 1])
std.err_ML
## [1] 0.02465517
p_value_ML <- 2 * pnorm(-abs(poly_ML$rho / std.err_ML))
p_value_ML
##              
## 1.927146e-94

并回答一个您没有问的重要问题：您会希望始终使用最大似然版本 ( ML=TRUE)，因为它更准确，除非您的计算机非常慢，在这种情况下，默认的两步方法是可以接受。