2016-04-26 9 views
0

は、私は図書館(MASS)パッケージからUScrimeデータの階層的回帰(のコードを実行するために、ぎざぎざで一緒にモデルを入れています。犯罪率が15の予測因子と応答、という。階層的回帰UScrime R /ぎざぎざ

私は私が解決することができません私のコードでいくつかのバグを持って

ぎざぎざモデル:。

model{ 
    for (i in 1:m){ 
     for (j in 1:47){ 
     y[j]~dnorm(beta[i,1]+beta[i,2]*x[j]+beta[i,3]*x[j]+beta[i,4]*x[j] 
     +beta[i,5]*x[j]+beta[i,6]*x[j]+beta[i,7]*x[j]+beta[i,8]*x[j] 
     +beta[i,9]*x[j]+beta[i,10]*x[j]+beta[i,11]*x[j]+beta[i,12]*x[j] 
     +beta[i,13]*x[j]+beta[i,14]*x[j]+beta[i,15]*x[j],invsig2)  
     } 
     beta[i,1:15]~dmnorm(theta[],invSig[,]) 
     } 
    theta[1:15]~dmnorm(mu0[],Lam.inv[,]) 
    invSig[1:15,1:15]~dwish(Lam[,],30) 
    invsig2~dgamma(.5,.5*sig20) 
    sig2<-1/invsig2 
    Sig<-inverse(invSig) 
} 

そしてRのコード:

crime=read.table("crime.dat",header=T) 
crime=as.matrix(crime) 

y=crime[,1];M=crime[,2];So=crime[,3];Ed=crime[,4];Po1=crime[,5];Po2=crime[,6] 
LF=crime[,7];M.F=crime[,8];Pop=crime[,9];NW=crime[,10];U1=crime[,11] 
U2=crime[,12];GDP=crime[,13];Ineq=crime[,14];Prob=crime[,15] 
Time=crime[,16] 
Crime=crime[,1] 
m <- length(unique(Crime)) # m = length of crimes 
n <- rep(NA,m) # setting up placeholder vectors 

for (j in 1:m){ 
n[j] <- sum(Crime==j) 
} 

beta_lse=matrix(nrow=m,ncol=15) 
sig2_lse=rep(NA,m) 

for (i in 1:m){ 
crime=y[1:47] 
x1=M[1:47];x2=So[1:47];x3=Ed[1:47];x4=Po1[1:47];x5=Po2[1:47];x6=LF[1:47] 
x7=M.F[1:47];x8=Pop[1:47];x9=NW[1:47];x10=U1[1:47];x11=U2[1:47] 
x12=GDP[1:47];x13=Ineq[1:47];x14=Prob[1:47];x15=Time[1:47] 
lse=lm(crime ~ 0+x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13+x14+x15) 
temp=anova(lse) 
beta_lse[i,]=c(lse[[1]]) 
sig2_lse[i]=temp[[3]][2]} 

Lam=cov(beta_lse) 
Lam.inv=solve(Lam) 
mu0=colMeans(beta_lse) 
sig20=mean(sig2_lse) 

library(rjags) 

forJags<-list(y=y,m=m,x=x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13+x14+x15, 
    Lam=Lam, Lam.inv=Lam.inv, mu0=mu0, sig20=sig20) 

inits<-list(beta=matrix(0,nrow=m,ncol=15),theta=c(0,0), invSig=.1*diag(2), invsig2=1) 

mod<-jags.model(file="CrimeHierarchicalRegression.bug",data=forJags,inits=inits) 
out<-coda.samples(model=mod,variable.names=c("theta","Sig","sig2","beta"),n.iter=10000,thin=5) 
summary(out) 

2つの太字のエラーError1とError2があり、それぞれの行でエラーが発生しています。最初のエラーは、「置き換えるアイテムの数が置換の長さの倍数ではない」と伝えます。 2番目のエラーは、その行のどこかに「予期しないシンボルがあります」と伝えています。

これらのエラーの原因が分かっている場合は、教えてください。前もって感謝します。

+0

私はこれが割り当てのためではなく、むしろ私たちの最終試験の準備として与えられた練習問題です。 – PapaSwankey

答えて

0

最初にエラーが発生した場合は、長さ15行のベクトルを16列の行列に入れようとしています。 2番目のエラーでは、x12=x12x13=x13の間にカンマはありません。

コンマを追加してbeta_lse=matrix(nrow=m,ncol=16)beta_lse=matrix(nrow=m,ncol=15)に変更してください。あなたはうまくいくはずです。もちろん、モデルの線形予測子で現在指定しているモデルにインターセプトを適合させたくないと仮定しています。その2つの特定のエラーを修正する必要があります。

+0

私はもはやそれらのエラーを取得していません。しかし、私が使っている行では、Lam.inv = solve(Lam)です。私は今、エラーを取得しています: "Lapackルーチンdgesv:システムは正確に単数形です:U [1,1] = 0"。私はそのコードを見て、それは基本的に私のラム行列が可逆ではないことを意味します。それがどのように是正されるかについてのアイデアはありますか? – PapaSwankey

+0

これは、あなたの行列が可逆ではなく、スタックオーバーフローの他の部分について議論されていることを意味します。 http://stackoverflow.com/questions/6572119/r-solvesystem-is-exactly-匿名 –

関連する問題