Simple variable intercept model help

Modeling
1

Hi,

I’m just trying to get my head around building models, and this is not working. I think you can see what I’m trying to do and I’m pretty sure the problem is in the loop in the transformed parameters block. Data setup with lm and plots up front. Thanks for your help:

library(tidyverse)
library(rstan)
library(shinystan)


n <- 50
j <- 3
beta <- 2
#noise <- rnorm(length(x), mean=0, sd=10)
x <- (1:n)
account <- (1:j)*20;
likeData <- data.frame(y = as.numeric(), x = as.numeric(), account = as.numeric())
for(i in 1:j){
   y = account[i] + beta*x + rnorm(length(x), mean=0, sd=10)
   tmpDF <- data.frame(y=y, x=x, account = account[i])
   likeData <- rbind(likeData, tmpDF)
}
likeData$account <- as.factor(likeData$account)
p <- ggplot(likeData, aes(x, y, color = account))+geom_point() 
p

lmMod <- lm(y ~ x + account, data=likeData)
summary(lmMod)

ggplot(likeData, aes(x=account, y=y))+geom_boxplot()

stanModel <- '
data {
  int<lower = 0> N; 
  int<lower = 0> J; 
  vector[N] account; 
  vector[N] y; 
  vector[N] x; 
} 
parameters { 
  vector[J] alpha; 
  real beta; 
  real alphaMu; 
  real alphaSigma;
  real ySigma;
} 
transformed parameters { 
  vector[N] yMu;
  for(i in 1:J) {
  yMu = alpha[i] + beta*x; 
  }
} 
model { 

  // prior

  alphaSigma ~ uniform(0,100);
  alpha ~ normal(alphaMu, alphaSigma);
  beta ~ normal(0, 10);

  // likelihood  

  ySigma ~ uniform(0,100);
  y ~ normal(yMu, ySigma); 

}
'

stanData <- list(
  N=nrow(likeData),
  J=nlevels(likeData$account),
  account = as.integer(likeData$account),
  y=likeData$y,
  x=likeData$x
)


fit <- stan(model_code = stanModel, data = stanData)

[Edit: Escaped code]

2

I’m afraid I odn’t see what you’re trying to do. Why do you keep reassigning yMu in the loop? It will wind up equal to alpha[J] + beta * x.

You also need to declare real<lower = 0, upper = 100> for alphaSigma and ySigma–support needs to match constraints. But we strongly discourage this prior for both statistical and computational reasons.

Also, you will need a non-centered parameterization for alpha. See the manual.

3

Thanks Mr. Carpenter,

It appears that lot of the postings out there about Stan originated before the program was vectorized? Going back to your manual, I fed in a matrix of predictors and it works fine (! once I dropped the intercept - I was feeding in a design matrix instead of a data matrix, just based on… oh, I don’t know, I just thought I needed to tell it there was an intercept to be considered for some reason). Sorry about taking up your time - RTFM :).

4

Yes, we’re still catching up on doc. And we need to clean up many of our example models. It’s not just vectorization—we’ve also added more convenient and stable compound distributions and functions. As well as changed a bit of the language to allow multi-indexing.