I'm trying to fit a categorical-logistic model using the painters dataset contained in the MASS library.
I divided the dataset in two parts, so i can predict in the future the values of School variable using the baseline category logistic multinomial model to assess the relationship between the categorical response School and the continuous explanatory variables. Unfortunately, my code runs or errors without even attempting a prediction.
library(MASS)
set.seed(43)
painters_new<-painters[sample(nrow(painters)) , ] #randomizing rows
painters_pred<-painters_new[50:54,]
painters_new<-painters_new[1:49,]
model<-function(){
#likelihood
for( i in 1:N){
y[i] ~ dcat(p[i,])
#formulating probabilities
p[i,1]<-ebx[i,1]/ebxsum[i]
p[i,2]<-ebx[i,2]/ebxsum[i]
p[i,3]<-ebx[i,3]/ebxsum[i]
p[i,4]<-ebx[i,4]/ebxsum[i]
p[i,5]<-ebx[i,5]/ebxsum[i]
p[i,6]<-ebx[i,6]/ebxsum[i]
p[i,7]<-ebx[i,7]/ebxsum[i]
p[i,8]<-ebx[i,8]/ebxsum[i]
#formulating denominator of probabilities
ebxsum[i]<-ebx[i,1]+ebx[i,2]+ebx[i,3]+ebx[i,4]+ebx[i,5]+ebx[i,6]+ebx[i,7]+ebx[i,8]
#formulating numerator of probabilities
ebx[i,1]<-exp(beta0[1]+beta1[1]*composition[i]+beta2[1]*drawing[i]+beta3[1]*colour[i]+beta4[1]*expression[i])
ebx[i,2]<-exp(beta0[2]+beta1[2]*composition[i]+beta2[2]*drawing[i]+beta3[2]*colour[i]+beta4[2]*expression[i])
ebx[i,3]<-exp(beta0[3]+beta1[3]*composition[i]+beta2[3]*drawing[i]+beta3[3]*colour[i]+beta4[3]*expression[i])
ebx[i,4]<-1 #baseline is Category 4
ebx[i,5]<-exp(beta0[5]+beta1[5]*composition[i]+beta2[5]*drawing[i]+beta3[5]*colour[i]+beta4[5]*expression[i])
ebx[i,6]<-exp(beta0[6]+beta1[6]*composition[i]+beta2[6]*drawing[i]+beta3[6]*colour[i]+beta4[6]*expression[i])
ebx[i,7]<-exp(beta0[7]+beta1[7]*composition[i]+beta2[7]*drawing[i]+beta3[7]*colour[i]+beta4[7]*expression[i])
ebx[i,8]<-exp(beta0[8]+beta1[8]*composition[i]+beta2[8]*drawing[i]+beta3[8]*colour[i]+beta4[8]*expression[i])
}
#priors
beta0[1:3]~dnorm(0,1.0E-06)
beta1[1:3]~dnorm(0,1.0E-06)
beta2[1:3]~dnorm(0,1.0E-06)
beta3[1:3]~dnorm(0,1.0E-06)
beta4[1:3]~dnorm(0,1.0E-06)
beta0[4]<-0
beta1[4]<-0
beta2[4]<-0
beta3[4]<-0
beta4[4]<-0
beta0[5:8]~dnorm(0,1.0E-06)
beta1[5:8]~dnorm(0,1.0E-06)
beta2[5:8]~dnorm(0,1.0E-06)
beta3[5:8]~dnorm(0,1.0E-06)
beta4[5:8]~dnorm(0,1.0E-06)
}
library(R2OpenBUGS)
model.file <- file.path(tempdir(),"model1.txt")
write.model(model, model.file)
y<-as.numeric(painters_new$School)
N<-49
composition<-painters$Composition
drawing<-painters$Drawing
colour<-painters$Colour
expression<-painters$Expression
data<-list("y","N","composition","drawing","colour","expression")
params<-c('beta0','beta1','beta2','beta3','beta4')
inits<-function(){list( beta0=rep(0,8),beta1=rep(0,8),beta2=rep(0,8),beta3=rep(0,8),beta4=rep(0,8))}
out<-bugs(data,inits,params,model.file,n.chains = 3
,n.iter=6000,codaPkg = TRUE,n.burnin = 1000)
Opening the log.txt file i get the following error message:
model is syntactically correct
data loaded
expected multivariate node
Obviously, this corresponds with a data structure error ' some variable that I use is expected to have different dimensions. Though I've tried different ways of defining the entities that I use, I'm completely sure that, according to theory, the nodes contained have the proper form. Is there any possibility that BUGS can't work with categorical and expects me to scale it to multinomial with 1 iteration?
Any ideas?
after various trial and modification, it seems that the problem's source was the way i'm trying to define beta-priors. BUGS does not support usual array tricks of R. I tried some if statements to dribble my way through and i figured that the usual way of an if statement in BUGS language(lol?) is using 2 functions of BUGS equals(x,z) and step(). So i did my homework about them and rearranged the whole prior
thing. model<-function(){
#likelihood
for( i in 1:N){
y[i] ~ dcat(p[i,])
#formulating probabilities
p[i,1]<-ebx[i,1]/ebxsum[i]
p[i,2]<-ebx[i,2]/ebxsum[i]
p[i,3]<-ebx[i,3]/ebxsum[i]
p[i,4]<-ebx[i,4]/ebxsum[i]
p[i,5]<-ebx[i,5]/ebxsum[i]
p[i,6]<-ebx[i,6]/ebxsum[i]
p[i,7]<-ebx[i,7]/ebxsum[i]
p[i,8]<-ebx[i,8]/ebxsum[i]
#formulating denominator of probabilities
ebxsum[i]<-ebx[i,1]+ebx[i,2]+ebx[i,3]+ebx[i,4]+ebx[i,5]+ebx[i,6]+ebx[i,7]+ebx[i,8]
#formulating numerator of probabilities
ebx[i,1]<-exp(beta0[1]+beta1[1]*composition[i]+beta2[1]*drawing[i]+beta3[1]*colour[i]+beta4[1]*expression[i])
ebx[i,2]<-exp(beta0[2]+beta1[2]*composition[i]+beta2[2]*drawing[i]+beta3[2]*colour[i]+beta4[2]*expression[i])
ebx[i,3]<-exp(beta0[3]+beta1[3]*composition[i]+beta2[3]*drawing[i]+beta3[3]*colour[i]+beta4[3]*expression[i])
ebx[i,4]<-1 #baseline is Category 4
ebx[i,5]<-exp(beta0[5]+beta1[5]*composition[i]+beta2[5]*drawing[i]+beta3[5]*colour[i]+beta4[5]*expression[i])
ebx[i,6]<-exp(beta0[6]+beta1[6]*composition[i]+beta2[6]*drawing[i]+beta3[6]*colour[i]+beta4[6]*expression[i])
ebx[i,7]<-exp(beta0[7]+beta1[7]*composition[i]+beta2[7]*drawing[i]+beta3[7]*colour[i]+beta4[7]*expression[i])
ebx[i,8]<-exp(beta0[8]+beta1[8]*composition[i]+beta2[8]*drawing[i]+beta3[8]*colour[i]+beta4[8]*expression[i])
}
#priors construction step-like
#f functions
for(j in 1:8){
f_branch[j,1]~dnorm(0,1.0E-6)
f_branch[j,2]<-0
#decision
if_branch[j]<-1+equals(j,4)
beta0[j]<-f_branch[j,if_branch[j]]
beta1[j]<-f_branch[j,if_branch[j]]
beta2[j]<-f_branch[j,if_branch[j]]
beta3[j]<-f_branch[j,if_branch[j]]
beta4[j]<-f_branch[j,if_branch[j]]
}}
library(R2OpenBUGS)
model.file <- file.path(tempdir(),"cat_log.txt")
write.model(model, model.file)
y<-as.numeric(painters_new$School)
N<-49
composition<-painters$Composition
drawing<-painters$Drawing
colour<-painters$Colour
expression<-painters$Expression
data<-list("y","N","composition","drawing","colour","expression")
params<-c('beta0','beta1','beta2','beta3','beta4')
out<-bugs(data,inits=NULL,params,model.file,n.chains = 3
,n.iter=6000,codaPkg = TRUE,n.burnin = 1000)
So now after dribbling the next set of errors, i get an error in R while log.txt file doesn't give me an indication of error type.
model is syntactically correct
data loaded
model compiled
Any ideas?
