######################################################
#Lecture 5 Script 
#Bayesian data analysis
#Anastasia Ejova 10/11/2014
#R101/PSY532
######################################################

#Useful resources
#Blog post re p-values in bayesglm output
  #http://stats.stackexchange.com/questions/31492/help-me-understand-p-values-in-bayesian-glm
#useful powerpoint:
  #http://rstudio-pubs-static.s3.amazonaws.com/15619_15a956c6bda84c3f97ef32793895689a.html#/10
#code for bayesglm and MCMCpack
  #http://www.unt.edu/rss/class/Jon/R_SC/Module10/BayesianGLM.R

#install.packages("arm")

library(arm)
library(ggplot2)

#ANOVA
bayesanova <- bayesglm(PostNaturalIoC ~ SeqCond * PreDBC_IOC, data = SS) #Cauchy prior (default)
  #See help file for how to set up a normal prior using prior.df
summary(bayesanova)

simulates <- coef(sim(bayesanova))
head(simulates)

posterior_Asc_vs_Desc <- simulates[,3]
hist(posterior_Asc_vs_Desc)
qplot(geom = "density", x = posterior_Asc_vs_Desc)
quantile(posterior_Asc_vs_Desc, c(.025, .975))

#(Intercept)                  0.25757    0.60240   0.428  0.66924 simulates[,1]   
#SeqCondU-shaped             -1.37424    0.85198  -1.613  0.10771 simulates[,2]
#SeqCondAscending            -1.78320    0.95976  -1.858  0.06408 . simulates [,3]
#SeqCondFlat                 -1.23566    0.85710  -1.442  0.15035   simulates[,4] 
#PreDBC_IOC                   0.06584    0.02066   3.186  0.00158 ** simulates[,5]
#SeqCondU-shaped:PreDBC_IOC   0.07016    0.02934   2.391  0.01736 *  simulates[,6]
#SeqCondAscending:PreDBC_IOC  0.09128    0.03253   2.806  0.00532 ** simulates[,7]
#SeqCondFlat:PreDBC_IOC       0.06173    0.02921   2.113  0.03535 *  simulates[,8]

#ANCOVA
bayesancova <- bayesglm(PostNaturalIoC ~ SeqCond + PreDBC_IOC, data = SS) 
summary(bayesancova)

simulatesanc <- coef(sim(bayesancova))
head(simulatesanc)

posterior_Asc_vs_Desc <- simulatesanc[,3]
hist(posterior_Asc_vs_Desc)
qplot(geom = "density", x = posterior_Asc_vs_Desc)
quantile(posterior_Asc_vs_Desc, c(.025, .975))

#Simulations in non-bayesian analysis (e.g., generalized linear modelling)
library(MASS)
ancova <- glm.nb(PostNaturalIoC ~ SeqCond + PreDBC_IOC, data = SS)
summary(ancova)
simulates0 <- coef(sim(ancova))
head(simulates0)
simulates_Asc_vs_Desc0 <- simulates0[,3]
hist(simulates_Asc_vs_Desc0)
qplot(geom = "density", x = posterior_Asc_vs_Desc0)
quantile(simulates_Asc_vs_Desc0, c(.025, .975))

#So what we are drawing is not the ACTUAL posterior displaying Markov Chain Monte Carlo (MCMC) 
#results. The density plots above were actually plots of regression slope simulations
#based on what is known as a MAP estimation procedure (see "Maximum a posteriori" estimation on
#Wikipedia). As the Environment window shows, simulates, simulatesanc and simulates0 contain only 
#100 values. An MCMC sample would contain thousands of cases (approximately 10 000).

#################################################
#Generalized linear modelling
#library(MASS)
glmnb <- glm.nb(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType + 
                  PreDBC_Total, data = SS) #From Lecture 4
summary(glmnb) #Note down value of dispersion parameter
#dispersion parameter(theta) is 2.4481 in this case

bayes.glmnb <- bayesglm(PostHowManySingleWins ~ SeqCond*PostHowManySingleCaptionType + 
                          PreDBC_Total, family = negative.binomial(theta = 2.4481), data = SS)
      #add value of dispersion parameter into the family argument
summary(bayes.glmnb)

simulatesglmnb <- coef(sim(bayes.glmnb)) #proceed with simulation as for ANOVA and ANCOVA

simulates_Asc_vs_Desc_glm <- simulatesglmnb[,3]
quantile(simulates_Asc_vs_Desc_glm, c(.025, .975))

simulates_Misses_vs_Goals_glm <- simulatesglmnb[,5]
quantile(simulates_Misses_vs_Goals_glm, c(.025, .975))

simulates_DescGoals_vs_FlatMisses_glm <- simulatesglmnb[,8]
quantile(simulates_DescGoals_vs_FlatMisses_glm, c(.025, .975))
#etc.

#################################################
#Logistic regression
#The next 10 lines are similar to Lecture 4
CG1$SumPreds <- as.numeric(CG1$PathCat) +
  as.numeric(CG1$PlayFreq) + as.numeric(CG1$Online)
CG2<- CG1[!is.na(CG1$SumPreds),]

summary(CG2$PlayFreq) #Categories not even, so will reduce number of categories.
#Will make a big category "Weekly or more" containing 53 people
library(plyr)
?mapvalues
CG2$PlayFreq <- mapvalues(CG2$PlayFreq, from = c("1-2 times per week", "3-4 times per week",
                                                     "Almost every day"), 
                            to = c("Weekly", "Weekly", "Weekly"))

bayeslogr <- bayesglm(PathCat ~ PlayFreq*Online, 
                      family=binomial(link="logit"), data = CG2)
        #again, Cauchy prior default
        #see this webpage and the help file for t-prior and normal prior options:
        #http://www.stat.columbia.edu/~gelman/standardize/bayesglm.R
display(bayeslogr)
summary(bayeslogr)
exp(bayeslogr$coefficients) #coefficients in odds ratio form

simulateslog <- coef(sim(bayeslogr))

post_Once_vs_LessThanMonthly <- simulateslog[,2] #Second line in summary(bayeslogr)
quantile(post_Once_vs_LessThanMonthly, c(.025, .975))

post_Once_vs_OncePerMonth <- simulateslog[,3]
quantile(post_Once_vs_OncePerMonth, c(.025, .975))

post_Once_vs_2to3PerMonth <- simulateslog[,4]
quantile(post_Once_vs_2to3PerMonth, c(.025, .975))

post_Once_vs_Weekly <- simulateslog[,5]
quantile(post_Once_vs_Weekly, c(.025, .975))

post_OnlineNo_vs_OnlineYes <- simulateslog[,6]
quantile(post_OnlineNo_vs_OnlineYes, c(.025, .975))

post_InteractionContrast1 <- simulateslog[,7]
quantile(post_InteractionContrast1, c(.025, .975))

post_InteractionContrast2 <- simulateslog[,8]
quantile(post_InteractionContrast2, c(.025, .975))

post_InteractionContrast3 <- simulateslog[,9]
quantile(post_InteractionContrast3, c(.025, .975))

post_InteractionContrast4 <- simulateslog[,10]
quantile(post_InteractionContrast4, c(.025, .975))

