我正在处理金融/经济数据,以防你想知道下面一些系数的大尺寸…我的一般问题与从R中的线性随机效应模型输出的参数系数的模拟有关.我是尝试使用来自R中相同模型的模型系数和方差 – 协方差(VCOV)矩阵生成β系数的随机样本.我的问题是：为什么我使用以下方法收到关于预期值的平方根的误差？来自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:
Warning message: In sqrt(ev$values) : NaNs produced
当我谷歌以下搜索字符串：“rmvnorm,”警告消息：在sqrt(ev $值)：NaNs产生,“我看到：
http://www.nickfieller.staff.shef.ac.uk/sheff-only/mvatasksols6-9.pdf在第4页上,此错误表示“负特征值”.虽然,我从概念上或实际上都不知道负特征值是什么,或者为什么在这种情况下会产生它们.
第二个搜索结果：[http://www.r-tutor.com/r-introduction/basic-data-types/complex2表示由于尝试取-1的平方根而产生此错误,这是“非复数值”(您不能取-1的平方根).
问题仍然存在,随着贝塔的随机生成,这里发生了什么,以及如何纠正这个问题？
sessionInfo() R version 3.0.2 (2013-09-25) Platform:
x86_64-apple-darwin10.8.0 (64-bit)Using the following packages/versions
mvtnorm_0.9-9994,
lme4_1.1-5,
Rcpp_0.10.3,
Matrix_1.1-2-2,
lattice_0.20-23
最佳答案 你的特征值有很大的尺度范围：
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
编辑：OP指出使用方法=“svd”与rmvnorm也有效.
如果您打印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))以完全精度保存方差 – 协方差矩阵更安全.
>如果你不知道什么是“肯定的”意味着你可能想要阅读它.关于covariance matrices和positive definite matrices的维基百科文章对你来说可能有点太技术性; this question on StackExchange更接近,但仍然有点技术性.我的Google之旅的下一个条目是this one,看起来是正确的.
