Hierarchical Regression UScrime R/JAGS

Question

I am putting together a model in JAGS to run code for a hierarchical regression of UScrime data (from the library(MASS) package. Crime rates being the response, with 15 predictors.

I have a few bugs in my code that I am unable to remedy.

JAGS Model:

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)
}

And the code in 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)

There are two bolded Errors, Error1 and Error2, with their respective lines causing the errors. The first error tells me that "number of items to replace is not a multiple of replacement length". The second error tells me that "there is an unexpected symbol" somewhere on that line.

Please let me know if you see how these errors are being caused. Thanks in advance.

I'd like to add that this is not for an assignment, but rather it is a practice problem given in preparation for our final exam. — PapaSwankey, Apr 26 '16 at 17:42

score 0 · Answer 1 · answered Apr 27 '16 at 14:21

0

For your first error you are trying to put vector of length 15 row wise into a matrix with 16 columns. For your second error, you do not have a comma between x12=x12 and x13=x13.

Add that comma and change beta_lse=matrix(nrow=m,ncol=16) to beta_lse=matrix(nrow=m,ncol=15) and you should be good to go. This is of course assuming that you do not want to fit an intercept to the model, which is what you have currently specified in the linear predictor of your model. That should fix those two specific errors.

answered Apr 27 '16 at 14:21

mfidino

3,030
1
9
13

I am no longer getting those errors. However, on the line where I am using Lam.inv = solve(Lam). I am now getting the error: " Lapack routine dgesv: system is exactly singular: U[1,1] = 0". I looked that code up, and it basically means that my Lam matrix is not invertible. Any ideas on how that is remedied? – PapaSwankey Apr 27 '16 at 16:23
This means that your matrix is not invertible, and has been discussed on other parts of stack overflow. http://stackoverflow.com/questions/6572119/r-solvesystem-is-exactly-singular – mfidino Apr 28 '16 at 13:29

Hierarchical Regression UScrime R/JAGS

1 Answers1