https://statmath.wu.ac.at/~gruen/BayesMix/bayesmix-intro.pdf

Using package bayesmix
##########################################################

library(bayesmix)

df <- read.table("D:/Dropbox/edps 590BAY/Lectures/7 Multiple Linear Regression/getting_what_paid_for_data.txt", header=TRUE)
names(df)
head(df)

sat <- as.data.frame(df$ave_tot)

model <- BMMmodel(sat, k = 2, initialValues = list(S0 = 2),
                  priors = list(kind = "independence",
                           hierarchical="tau"))
   
control <- JAGScontrol(variables = c("mu", "tau", "eta", "S"),
 burn.in = 1100, n.iter = 40000, seed = 10)
 
z <- JAGSrun(sat, model = model, control = control)

plot(z)

print(z)

#classification

results <- z$results

class <- matrix(NA ,nrow=50, ncol=2)
for (state in 1:50) {
   c <- table(results[,state])
   class[state,] <- c
}

row.names(class) <- df$state

df$postprob_SAT <- class[,1]/50000
df$postprob_ACT <- class[,2]/50000

plot(df$ave_tot, df$postprob_SAT , pch="n", ylab="Probability class 1",
         main="Class 1 (SAT)")
text(df$ave_tot, df$postprob_c1, labels=df$state)

par(mfrow=c(1,1))
plot(df$postprob_SAT, df$ave_tot , type="p", 
       ylab="Average total SAT",
       xlab="Posterior Probability SAT",
       main="Classification of states")
text(df$postprob_SAT, df$ave_tot , df$state_abv)


d <- df[, c(1,8,11,12)]

(dorder <- d[order(d$postprob_SAT),])