Bayes.regression = function(y,X,b0,B0,nu0,s20){
  p      = ncol(X)
  b0     = rep(b0,p)
  B0     = diag(B0,p)
  nu0s20 = nu0*s20
  iB0    = solve(B0)
  iB0b0  = iB0%*%b0
  B1     = solve(iB0+t(X)%*%X)
  tcB1   = t(chol(B1))
  b1     = B1%*%(iB0b0+t(X)%*%y)
  nu1    = nu0 + (n-p)
  s21    = (nu0s20 + t(y-X%*%b1)%*%y + t(b0-b1)%*%iB0b0)/nu1
  se     = sqrt(diag(s21[1,1]*B1))
  m      = X%*%b0
  V      = s20*(diag(1,n)+X%*%B0%*%t(X))
  res    = t(chol(solve(V)))%*%(y-m)
  pred   = sum(dt(res,df=nu0,log=TRUE))
  return(list(b1=b1,B1=B1,nu1=nu1,s21=s21,se=se,pred=pred))
}


data = read.table("eleicaopresidencial-2022-mg.txt",header=TRUE)

attach(data)

n  = nrow(data)
y  = log(votoslula/(1-votoslula))
x1 = abstencoes
x2 = log10(populacao)

par(mfrow=c(1,3))
plot(x1,y,xlab="Abstenções",ylab="Log odds para Lula")
title(paste("correlação=",round(cor(x1,y),2),sep=""))
plot(x2,y,xlab="População (log base 10)",ylab="Log odds para Lula")
title(paste("correlação=",round(cor(x2,y),2),sep=""))
plot(x1,x2,xlab="Abstenções",ylab="População (log base 10)")
title(paste("correlação=",round(cor(x1,x2),2),sep=""))

#################################
# Classical linear regression
#################################
fit.M0 = lm(y~x1)
fit.M1 = lm(y~x2)
fit.M2 = lm(y~x1+x2)

summary(fit.M0)
summary(fit.M1)
summary(fit.M2)

r2adj = c(
summary(fit.M0)$adj.r.squared,
summary(fit.M1)$adj.r.squared,
summary(fit.M2)$adj.r.squared)

# Deriving yourself
X        = cbind(1,x1,x2)
p        = ncol(X)
beta.mle = solve(t(X)%*%X)%*%t(X)%*%y
s.mle    = sqrt(sum((y-X%*%beta.mle)^2)/(n-p))
se.mle   = s.mle*sqrt(diag(solve(t(X)%*%X)))

cbind(beta.mle,se.mle)
s.mle



#################################
# Bayesian linear regression
#################################
b0     = 0
B0     = 4
nu0    = 10
s20    = 1


X1 = cbind(1,x1)
X2 = cbind(1,x2)
X3 = cbind(1,x1,x2)

fit.bayes1 = Bayes.regression(y,X1,b0,B0,nu0,s20)
fit.bayes2 = Bayes.regression(y,X2,b0,B0,nu0,s20)
fit.bayes3 = Bayes.regression(y,X3,b0,B0,nu0,s20)

fit.bayes1$pred
fit.bayes2$pred
fit.bayes3$pred

# Comparing MLE and posterior analysis
cbind(beta.mle,se.mle)
cbind(fit.bayes3$b1,fit.bayes3$se)

# Model comparison:  R2, predictive
# Bayes factor and posterior model probability
pred = c(fit.bayes1$pred,fit.bayes2$pred,fit.bayes3$pred)
pmp  = exp(pred-max(pred))
pmp  = pmp/sum(pmp)
cbind(r2adj,pred,pmp)