1

我需要分析一些具有以下结构的模拟数据：

h   c   x1              y1              x1c10
1   0   37.607056431    104.83097593    5
1   1   27.615251557    140.85532974    10
1   0   34.68915314     114.59312842    2
1   1   30.090387454    131.60485642    9
1   1   39.274429397    106.76042522    10
1   0   33.839385007    122.73681319    2
...

其中 h 的范围从 1 到 2500，并索引 Monte Carlo 样本，每个样本有 1000 个观测值。我正在使用以下代码分析这些数据，该代码为我提供了两个对象（fnN1、fdQB101）：

mc<-2500 ##create loop index
fdN1<-matrix(0,mc,1000)
fnQB101 <- matrix(0,mc,1000) ##create 2500x1000 storage matrices, elements zero

for(j in 1:mc){

fdN1[j,] <- dnorm(residuals(lm(x1 ~ c,data=s[s$h==j,])), 
                     mean(residuals(lm(x1 ~ c,data=s[s$h==j,]))), 
                          sd(residuals(lm(x1 ~ c,data=s[s$h==j,]))))

x1c10<-as.matrix(subset(s,s$h==j,select=x1c10))

fdQB100 <- as.matrix(predict(polr(as.factor(x1c10) ~ c , 
                                    method="logistic", data=s[s$h==j,]),
                                         type="probs"))

indx10<- as.matrix(cbind(as.vector(seq(1:nrow(fdQB100))),x1c10))

fdQB101[j,] <- fdQB100[indx10]

}

对象 fdN1 和 fdQB101 是 2500x1000 矩阵，以预测概率作为元素。我需要从这个循环中创建一个可以用 lapply() 或 mclapply() 调用的函数。当我将其包装在以下函数命令中时：

ndMC <- function(mc){

for(j in 1:mc){
...
}
return(list(fdN1,fdQB101))

}
lapply(mc,ndMC)

对象 fdN1 和 fdQB101 每个都作为 2500x1000 的零矩阵返回，而不是预测的概率。我究竟做错了什么？

4

1 回答 1

1

你应该可以用这个data.table包来做到这一点。这是一个例子：

library(data.table)
dt<-data.table(h=rep(1L,6), c=c(0L,1L,0L,1L,1L,0L),
           X1=c(37.607056431,27.615251557,34.68915314,30.090387454,39.274429397,33.839385007),
           y1=c(104.83097593,140.85532974,114.59312842,131.60485642,106.76042522,122.73681319),
           x1c10=c(5L,10L,2L,9L,10L,2L))

## Create a linear model for every grouping of variable h:
fdN1.partial<-dt[,list(lm=list(lm(X1~c))),by="h"]

## Retrieve the linear model for h==1:
fdN1.partial[h==1,lm]
## [[1]]
## 
## Call:
## lm(formula = X1 ~ c)
## 
## Coefficients:
## (Intercept)            c  
##      35.379       -3.052

您还可以编写一个函数来概括此解决方案：

f.dnorm<-function(y,x) {
  f<-lm(y ~ x)
  out<-list(dnorm(residuals(f), mean(residuals(f)), sd(residuals(f))))
  return(out)
}

## Generate two dnorm lists for every grouping of variable h:
dt.lm<-dt[,list(dnormX11=list(f.dnorm(X1,rep(1,length(X1)))), dnormX1c=list(f.dnorm(X1,c))),by="h"]

## Retrieve one of the dnorm lists for h==1:
unlist(dt.lm[h==1,dnormX11])
##          1          2          3          4          5          6 
## 0.06296194 0.03327407 0.08884549 0.06286739 0.04248756 0.09045784
于 2013-03-11T04:23:03.220 回答