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で、モンテカルロサンプルのインデックスを作成します。各サンプルには1000個の観測値があります。これらのデータを、2つのオブジェクト（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 に答える