R中的下标越界错误答案

【问题标题】：Subscript out of bound error in RR中的下标越界错误
【发布时间】：2015-04-03 23:03:16
【问题描述】：

使用 stats 包中的 factanal 函数进行因子分析时。

我尝试了以下事情。

library(mirt)
library(ltm)
library(psych)
library(stats)
data(SAT12)
data=SAT12
cor_mat=polychoric(data,  ML=TRUE, global=F)
 fit <- factanal(factors=2, n.obs=nrow(data), covmat=cor_mat$rho)


Divide_item_Factor_Loading(fit)

当我尝试运行 Divide_item_Factor_Loading(fit) 时出现一个名为

的错误

   Error in a[[i]][[2]] : subscript out of bounds

弹出。

我的 Divide_item_Factor_Loading 完整代码：

Divide_item_Factor_Loading=function(fit)
{
  a=list()
  items=NULL
  for(i in 1:nrow(fit$loadings)) ######corresponding to rows of loading matrix
  {
    k=which(fit$loadings[i,]==max(abs(fit$loadings[i,])))  
    a[[i]]=c(i,as.numeric(k))
  } 
  fact_item_mat=matrix(, nrow=nrow(fit$loadings), ncol=ncol(fit$loadings))
  for(j in 1:(ncol(fit$loadings)))
  {
    for(i in 1:(nrow(fit$loadings)))
    {
      if(a[[i]][[2]]==j) {fact_item_mat[i,j]=a[[i]][[1]]}
    }    
  }
  nam=names(fit$loadings[,1])
  factor=list()
  for(i in 1:ncol(fit$loadings))
  {
    factor[[i]]=sort(fact_item_mat[,i], decreasing = FALSE, na.last = NA)
    fac=factor[[i]]
    fac=nam[fac]
    factor[[i]]=fac
  }
  names(factor)=paste("factor", 1:ncol(fit$loadings), sep="")
  return(factor)
}

我现在应该采取什么步骤来避免这个错误？

【问题讨论】：

正如 Andrei 所说，这只是打印方法——所有值都在负载矩阵中。查看loadings(fit) 和loadings(fit)[]。由于打印方法，它们看起来不同，但加载矩阵是相同的。检查str(loadings(fit)) 和str(loadings(fit)[])。再次检查loadings(fit)[] == loadings(fit)
@Irri，在LSAT 数据上传递您的代码不会给我任何错误（尽管“ML 选项已从 psych 包中的 polychoric 函数中删除。请修复通话”）并做Divide_item_Factor_Loading(fit)，我得到$factor1 [1] "Item 5" $factor2 [1] "Item 1" "Item 2" "Item 3" "Item 4"。所以，我很抱歉，但我真的不明白你的问题......
@CathG，能否请您提供您的电子邮件 ID。这样我就可以给你有问题的数据集。
@Irri，我不希望它成为解决方案。你真的应该重新表述你的问题，因为你给出了 2 个有效的例子......相反，你应该给出一个最小的可重复的例子，导致你得到错误（并给出预期的输出）。 Here 是一些可以帮助你的建议。
@CathG，我再次编辑了我的问题。我希望这个可以帮助您重现我遇到的错误。

标签： r factor-analysis

【解决方案1】：

运行您的代码并调试您的函数（使用debug 函数）我知道您为什么遇到“下标越界”错误：

变量 a 的第 15 个元素（以及其他元素）的长度为 1，因此当您尝试访问 a[[15]][2] 时，R 不高兴...
这个元素的长度只有 1 而不是 2 的原因是因为因子的最大绝对值达到了负值，而您要问哪个值（非绝对值）等于这个最大绝对值，所以答案是否定的……

因此您需要将行
which(fit$loadings[i,]==max(abs(fit$loadings[i,]))) 更改为which(abs(fit$loadings[i,])==max(abs(fit$loadings[i,])))
你会得到：

Divide_item_Factor_Loading(fit)
#$factor1
 #[1] "Item.1"  "Item.4"  "Item.6"  "Item.7"  "Item.8"  "Item.9"  "Item.10" "Item.11" "Item.13" "Item.14" "Item.15"
#[12] "Item.17" "Item.19" "Item.20" "Item.24" "Item.26" "Item.27" "Item.28" "Item.29"

#$factor2
 #[1] "Item.2"  "Item.3"  "Item.5"  "Item.12" "Item.16" "Item.18" "Item.21" "Item.22" "Item.23" "Item.25" "Item.30"
#[12] "Item.31" "Item.32"

即使被调试的函数现在可以工作，我认为你应该改变它，因为这比它应该的更复杂。

我对替代功能的建议：

Divide_item_Factor_Loading_v2<-function(fit){
     a<-apply(fit$loadings,1,function(facs) which(abs(facs)==max(abs(facs))))
     return(list(factor1=names(a)[a==1],factor2=names(a)[a==2]))
}

这会为您的 fit 对象提供与您的（调试的）函数完全相同的结果：

Divide_item_Factor_Loading_v2(fit)
#$factor1
 #[1] "Item.1"  "Item.4"  "Item.6"  "Item.7"  "Item.8"  "Item.9"  "Item.10" "Item.11" "Item.13" "Item.14" "Item.15"
#[12] "Item.17" "Item.19" "Item.20" "Item.24" "Item.26" "Item.27" "Item.28" "Item.29"

#$factor2
 #[1] "Item.2"  "Item.3"  "Item.5"  "Item.12" "Item.16" "Item.18" "Item.21" "Item.22" "Item.23" "Item.25" "Item.30"
#[12] "Item.31" "Item.32"

【讨论】：

非常感谢！ :)
@Irri，不客气，记住debug 函数，当函数产生错误时它非常有用；-)

【解决方案2】：

要更改载荷的打印方式，请使用 cutoff 参数到 print.loadings。

试试这样的：

print(fit$loadings, cutoff=0)

实际矩阵包含所有值。

print(loadings(fit), cutoff=0)

Loadings:
       Factor1 Factor2
Item 1 0.014   0.418  
Item 2 0.130   0.350  
Item 3 0.036   0.553  
Item 4 0.166   0.294  
Item 5 0.990   0.125  

               Factor1 Factor2
SS loadings      1.025   0.705
Proportion Var   0.205   0.141
Cumulative Var   0.205   0.346

现在提取每个因子的最大负载，使用apply()：

apply(loadings(fit), 2, max)

  Factor1   Factor2 
0.9895743 0.5531770

【讨论】：

感谢您的回复。我不想打印它。我需要将该负载矩阵存储在一个对象中，并使用它来计算一个项目的最大负载量，正如我在编辑问题的上面代码中所示的那样。
好了，所有的载荷都被存储了。他们只是在打印时被压制。我将编辑我的答案。
感谢您的回复安德烈。但是我需要为每个项目存储具有最大负载的因子。这就是列表“a”正在做的事情。不只是每个因素的一个最大值。但不幸的是，它无法比较这些值（而且我认为这些值是不可见的，这就是它无法比较的原因）。我将再次编辑我的问题。请看一下。

【解决方案3】：

检查?loadings，您会发现有一个cutoff 参数定义了一个“小于此（绝对值）的负载被抑制”的值。

【讨论】：