Bayes

UCLAのInstitute for Digital Research and EducationにあるRによる順序ロジスティック回帰分析の例

https://stats.idre.ucla.edu/r/dae/ordinal-logistic-regression/

require(foreign)

require(MASS)

d <- read.dta("https://stats.idre.ucla.edu/stat/data/ologit.dta")

head(dat)

require(foreign)

require(MASS)

d <- read.dta("https://stats.idre.ucla.edu/stat/data/ologit.dta")

head(d)

result <- polr(apply ~ pared + public + gpa, data = d, Hess=TRUE)

summary(result)

# 通常の順序ロジスティック回帰分析

library(rstan)

ordered <- "

data {

int<lower=1> N; /// sample size

int<lower=1> K; /// Number of categories of the dependent variable

int<lower=1> M; /// Number of independent variables

matrix[N,M] x; /// independent variables

int<lower=1,upper=K> y[N]; /// dependent variable

}

parameters {

vector[M] beta;

ordered[K-1] c;

}

model {

vector[K] theta;

for (i in 1:N) {

theta[1] = 1-inv_logit(x[i]*beta - c[1]);

for (k in 2:(K-1))

theta[k] = inv_logit(x[i]*beta- c[k-1]) - inv_logit(x[i]*beta - c[k]);

theta[K] = inv_logit(x[i]*beta - c[K-1]);

y[i] ~ categorical(theta);

beta ~ normal(0,100);

theta ~ normal(0,100);

}

# Independent Vairables

x <- model.matrix(d$apply~d$pared + d$public + d$gpa-1,d)

fit <- sampling(model, data = list(N = length(d$apply), K = max(as.numeric(d$apply)), M=ncol(x),y = as.numeric(d$apply), x = x))

print(fit)

names(fit)

stan_trace(fit,par="beta")

stan_trace(fit,par="c")

stan_hist(fit,pars="beta")

stan_dens(fit,pars="beta",separate_chains = TRUE)

stan_ac(fit,pars="beta",separate_chains = TRUE)

# 一般化順序ロジット

library(rstan)

gen_ordered <- "

data {

int<lower=1> N; /// sample size

int<lower=1> J; /// Number of categories of the dependent variable

int<lower=1> K; /// Number of independent variables

matrix[N,K] x; /// independent variables data (proportional)

int<lower=1> L; /// Number of independent variables

matrix[N,L] z; /// independent variables data (not proportional)

int<lower=1,upper=J> y[N]; /// dependent variable

}

parameters {

vector[K] beta;

matrix[L,J-1] betaj;

ordered[J-1] c;

}

model {

vector[J] theta;

for (i in 1:N) {

theta[1] = 1 - inv_logit(x[i]*beta + z[i]*betaj[,1] - c[1]);

for (j in 2:(J-1))

theta[j] = inv_logit(x[i]*beta + z[i]*betaj[,j-1] - c[j-1]) - inv_logit(x[i]*beta + z[i]*betaj[,j] - c[j]);

theta[J] = inv_logit(x[i]*beta + z[i]*betaj[,J-1] - c[J-1]);

y[i] ~ categorical(theta);

theta ~ normal(0,100);

}

model <- stan_model(model_code = gen_ordered)

x <- model.matrix(d$apply ~ d$pared + d$public-1,d)

z <- model.matrix(d$apply ~ d$gpa-1,d)

fit <- sampling(model, data = list(N = length(d$apply),

J = max(as.numeric(d$apply)),K=ncol(x),L=ncol(z),y=as.numeric(d$apply),x=x,z=z),

chains=2,iter=400,warmup=200)