分层回归UScrime R/JAGS

我组建了一个模型中JAGS为UScrime数据的层次回归（从库（MASS）包运行代码。犯罪率作为响应，用15个预测因子。分层回归UScrime R/JAGS

我在我的代码中的一些错误，我无法弥补

JAGS型号：

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)

有两个加粗的错误，Error1和Error2，它们各自的行导致错误。第一个错误告诉我“要替换的项目数量不是替换长度的倍数”。第二个错误告诉我，“有一个意外的符号”在该行的某处。

请让我知道，如果你看到这些错误是如何造成的。提前致谢。

来源

2016-04-26 PapaSwankey

我想补充一点，这不是一项任务，而是在准备期末考试时给出的练习题。 – PapaSwankey

对于您的第一个错误，您试图将长度为15的行的矢量放入16列的矩阵中。对于第二个错误，在x12=x12和x13=x13之间没有逗号。

添加该逗号并将beta_lse=matrix(nrow=m,ncol=16)更改为beta_lse=matrix(nrow=m,ncol=15)并且您应该很好。这当然假设您不想为模型拟合截距，这是您目前在模型的线性预测器中指定的截距。这应该解决这两个具体的错误。

来源

2016-04-27 14:21:18

我不再收到这些错误。但是，在我使用Lam.inv = solve（Lam）的线上。我现在得到的错误是：“Lapack routine dgesv：system is exactly singular：U [1,1] = 0”。我查看了代码，这基本上意味着我的Lam矩阵是不可逆的。任何想法如何解决？ – PapaSwankey

这意味着您的矩阵是不可逆的，并且已经在堆栈溢出的其他部分进行了讨论。 http://stackoverflow.com/questions/6572119/r-solvesystem-is-exactly-singular –

分层回归UScrime R/JAGS

回答

相关问题