为什么rmvnorm()函数返回"In sqrt(ev$values) : NaNs produced"，这是什么错误，如何纠正或避免？

Question

为什么rmvnorm()函数返回"In sqrt(ev$values) : NaNs produced"，这是什么错误，如何纠正或避免？

4

我正在处理与金融/经济数据相关的内容，如果您想知道下面某些系数的巨大大小...我的一般问题涉及在R中模拟线性随机效应模型的参数系数输出。我试图使用模型系数和来自同一模型的方差-协方差（VCOV）矩阵生成beta系数的随机样本。我的问题是：为什么我收到有关使用mvtnorm{}包中的rmvnorm()函数的预期值平方根的错误的下面警告？我该如何处理这个警告/问题？

#Example call: lmer model with random effects by YEAR
#mlm<-lmer(DV~V1+V2+V3+V2*V3+V4+V5+V6+V7+V8+V9+V10+V11+(1|YEAR), data=dat)
#Note: 5 years (5 random effects total)

#LMER call yields the following information:
coef<-as.matrix(c(-28037800,0.8368619,2816347,8681918,-414002.6,371010.7,-26580.84,80.17909,271.417,-239.1172,3.463785,-828326))

sigma<-as.matrix(rbind(c(1834279134971.21,-415.95,-114036304870.57,-162630699769.14,-23984428143.44,-94539802675.96,
                       -4666823087.67,-93751.98,1735816.34,-1592542.75,3618.67,14526547722.87),
                 c(-415.95,0.00,41.69,94.17,-8.94,-22.11,-0.55,0.00,0.00,0.00,0.00,-7.97),
                 c(-114036304870.57,41.69,12186704885.94,12656728536.44,-227877587.40,-2267464778.61,
                       -4318868.82,8909.65,-355608.46,338303.72,-321.78,-1393244913.64),
                 c(-162630699769.14,94.17,12656728536.44,33599776473.37,542843422.84,4678344700.91,-27441015.29,
                       12106.86,-225140.89,246828.39,-593.79,-2445378925.66),
                 c(-23984428143.44,-8.94,-227877587.40,542843422.84,32114305557.09,-624207176.98,-23072090.09,
                       2051.16,51800.37,-49815.41,-163.76,2452174.23),
                 c(-94539802675.96,-22.11,-2267464778.61,4678344700.91,-624207176.98,603769409172.72,90275299.55,
                       9267.90,208538.76,-209180.69,-304.18,-7519167.05),
                 c(-4666823087.67,-0.55,-4318868.82,-27441015.29,-23072090.09,90275299.55,82486186.42,-100.73,
                       15112.56,-15119.40,-1.34,-2476672.62),
                 c(-93751.98,0.00,8909.65,12106.86,2051.16,9267.90,-100.73,2.54,8.73,-10.15,-0.01,-1507.62),
                 c(1735816.34,0.00,-355608.46,-225140.89,51800.37,208538.76,15112.56,8.73,527.85,-535.53,-0.01,21968.29),
                 c(-1592542.75,0.00,338303.72,246828.39,-49815.41,-209180.69,-15119.40,-10.15,-535.53,545.26,0.01,-23262.72),
                 c(3618.67,0.00,-321.78,-593.79,-163.76,-304.18,-1.34,-0.01,-0.01,0.01,0.01,42.90),
                 c(14526547722.87,-7.97,-1393244913.64,-2445378925.66,2452174.23,-7519167.05,-2476672.62,-1507.62,21968.29,
                        -23262.72,42.90,229188496.83)))
#Error begins here:
betas<-rmvnorm(n=1000, mean=coef, sigma=sigma)
#rmvnorm breaks, Error returned:

警告信息：在sqrt（ev $ values）中：产生NaNs

当我谷歌搜索以下搜索字符串：“rmvnorm，“警告消息：在sqrt（ev $ values）中：产生NaNs”时，我发现： http://www.nickfieller.staff.shef.ac.uk/sheff-only/mvatasksols6-9.pdf 第4页指出此错误表示“负特征值”。尽管我不知道什么是负特征值或为什么会在这种情况下产生。

第二个搜索结果：[http://www.r-tutor.com/r-introduction/basic-data-types/complex 2表明，此错误是由于试图取-1的平方根而引起的，这是“不是复杂值”（无法取-1的平方根）。

问题仍然存在，即随机生成beta值的情况是怎样的，如何进行更正？

sessionInfo() R版本3.0.2（2013-09-25）平台： x86_64-apple-darwin10.8.0（64位）

使用以下软件包/版本 mvtnorm_0.9-9994， lme4_1.1-5， Rcpp_0.10.3， Matrix_1.1-2-2， lattice_0.20-23

- DV Hughes

@Fernando：这些是vcov()函数返回的值。您建议如何确定“vcov矩阵值是否错误”？请详细说明？ - DV Hughes

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

您的特征值范围非常广:

range(eigen(sigma)$values)
## [1] -1.005407e-05  1.863477e+12

我更喜欢使用MASS包中的mvrnorm函数，因为它会随着R的安装自动安装。此外，它似乎更加健壮：

set.seed(1001)
m <- MASS::mvrnorm(n=1000, mu=coef, Sigma=sigma)  ## works fine

编辑：楼主指出使用rmvnorm和method="svd"也可以。

如果打印MASS::mvrnorm的代码，或者debug(MASS:mvrnorm)并逐步执行它，你会发现它使用了

if (!all(ev >= -tol * abs(ev[1L]))) stop("'Sigma' is not positive definite")

（其中ev是特征值向量，按降序排序，因此ev [1]是最大的特征值），以判断方差-协方差矩阵的正定性。在这种情况下，ev [1L]约为2e12，tol为1e-6，因此这将允许负特征值的数量达到约2e6。在这种情况下，最小特征值为-1e-5，远低于公差。

MASS :: mvrnorm更低，使用pmax（ev，0）--也就是说，如果它已经确定特征值未达到公差以下（即未通过上述测试），则只会将负值截断为零，对于实际应用来说应该是可以接受的。

如果您坚持使用rmvnorm，则可以使用Matrix :: nearPD，该函数尝试强制使矩阵为正定矩阵--它返回一个列表，其中包含（除其他内容外）特征值和“正定矩阵”：

m <- Matrix::nearPD(sigma)
range(m$eigenvalues)
## [1] 1.863477e+04 1.863477e+12

从矩阵中计算出的特征值并不完全相同——nearPD和eigen使用略有不同的算法——但它们非常接近。

range(eigen(m$mat)$values)
## [1] 1.861280e+04 1.863477e+12

更一般地说，

特征值范围巨大的部分原因可能是预测变量的比例尺度非常不同。如果可能的话，缩放输入数据以使方差更加相似（即使所有数值计算更加稳定）可能是一个好主意--在生成值后，您总是可以重新缩放值。
当矩阵非常接近奇异时（即某些特征值非常接近零），小的数值差异会改变特征值的符号。特别是，如果您复制并粘贴值，则可能会失去一些精度并导致此问题。使用dput（vcov（fit））或save（vcov（fit））以完全精度保存方差协方差矩阵更安全。
如果您不知道“正定”是什么意思，您可能需要了解一下。维基百科关于协方差矩阵和正定矩阵的文章可能有点太技术化了；这个StackExchange上的问题更接近，但仍然有点技术化。我在谷歌上的下一个链接是这个，看起来还不错。