RStan でカテゴリカル変数の単回帰

係数の和＝０の制約をいれる場合。

##Fit with categorical variable 2
x0<-sample(1:4,100,replace=T)
library(dummies)
x<-dummy(x0)
x1<-x[,1]
x2<-x[,2]
x3<-x[,3]
x4<-x[,4]

a1<-2
a2<--2
a3<-10
a4<--10

y<-a1*x1+a2*x2+a3*x3+a4*x4+rnorm(100,sd=0.03)

#library("rstan")

data<-list(y=y,x2=x2,x1=x1,x3=x3,x4=x4,N=100,K=4)
fit<-stan(file="lm4.stan",data=data,seed=1234)
summary(fit)$summary

Rstan lm4.stan

data{
  
  int N;
  int K;
  real y[N];
  real x1[N];
  real x2[N];
  real x3[N];
  real x4[N];

}


parameters {
  
  real a[K-1];
  real b;
  real <lower=0> sigma;
  real <lower=0,upper=0.1> sigma_a;
}

transformed parameters{
  real y_base[N];
  real aa[K];
  real sum1=0;


  for(i in 1:(K-1)){
    aa[i]=a[i];
    sum1=sum1+aa[i];
  }
  aa[K]=0-sum1;

  for(i in 1:N){
   
   
    y_base[i]=aa[1]*x1[i]+aa[2]*x2[i]+aa[3]*x3[i]+aa[4]*x4[i]+b;
    
    
  }
  
  
}

model {
  
  for(i in 1:N){
    y[i] ~ normal(y_base[i],sigma);
  }  
}

カテゴリが0の制約を入れる場合
a1x1+a2x2+a3x4+a4x4=0;

##Fit with categorical variable 2
x0<-sample(1:4,100,replace=T)
library(dummies)
x<-dummy(x0)
x1<-x[,1]
x2<-x[,2]
x3<-x[,3]
x4<-x[,4]

a1<-2
a2<--2
a3<-10
a4<--10

y<-a1*x1+a2*x2+a3*x3+a4*x4+rnorm(100,sd=0.03)


data<-list(y=y,x2=x2,x1=x1,x3=x3,x4=x4,N=100,K=4)
fit<-stan(file="lm5.stan",data=data,seed=1234)
summary(fit)$summary

ms<- rstan::extract(fit)

apply(ms$aa2,2,mean)

lm5.stan

data{
  
  int N;
  int K;
  real y[N];
  real x1[N];
  real x2[N];
  real x3[N];
  real x4[N];

}


parameters {
  
  real a[K-1];
  real b;
  real <lower=0> sigma;
  real <lower=0,upper=0.1> sigma_a;
}

transformed parameters{
  real y_base[N];
  real aa[K];
  real sum1=0;


  for(i in 1:(K-1)){
    aa[i]=a[i];
  //  sum1=sum1+aa[i];
  }
  aa[K]= -(aa[1]*sum(x1)+aa[2]*sum(x2)+aa[3]*sum(x3))/sum(x4);

  for(i in 1:N){
   
   
    y_base[i]=aa[1]*x1[i]+aa[2]*x2[i]+aa[3]*x3[i]+aa[4]*x4[i]+b;
    
    
  }
  
  
}

model {
  
  for(i in 1:N){
    y[i] ~ normal(y_base[i],sigma);
  }  
}