Fixed bug in dt

misc_code/jags_linear_regression_R_squared.R

@@ -0,0 +1,79 @@
+# Header ------------------------------------------------------------------
+
+# Fitting a linear regression in JAGS and calculate the R-squared value
+# Andrew Parnell
+
+# In this code we generate some data from a simple linear regression model and fit is using jags. We then compute an R-squared values
+
+# Some boiler plate code to clear the workspace, and load in required packages
+rm(list=ls()) # Clear the workspace
+library(R2jags)
+
+# Maths -------------------------------------------------------------------
+
+# Description of the Bayesian model fitted in this file
+# Notation:
+# y_i = repsonse variable for observation t=i,..,N
+# x_i = explanatory variable for obs i
+# alpha, beta = intercept and slope parameters to be estimated
+# sigma = residual standard deviation
+
+# Likelihood:
+# y[i] ~ N(alpha + beta * x[i], sigma^2)
+# Prior
+# alpha ~ N(0,100) - vague priors
+# beta ~ N(0,100)
+# sigma ~ U(0,10)
+
+# Simulate data -----------------------------------------------------------
+
+# Some R code to simulate data from the above model
+n = 100
+alpha = 2
+beta = 3
+sigma = 1
+# Set the seed so this is repeatable
+set.seed(123)
+x = sort(runif(n, 0, 10)) # Sort as it makes the plotted lines neater
+y = rnorm(n, mean = alpha + beta * x, sd = sigma)
+
+# Also creat a plot
+plot(x, y)
+lines(x, alpha + beta * x)
+
+# Jags code ---------------------------------------------------------------
+
+# Jags code to fit the model to the simulated data
+
+model_code = '
+model
+{
+  # Likelihood
+  for (i in 1:n) {
+    y[i] ~ dnorm(fits[i], sigma^-2)
+    fits[i] = alpha + beta * x[i]
+  }
+
+  # Priors
+  alpha ~ dnorm(0, 100^-2)
+  beta ~ dnorm(0, 100^-2)
+  sigma ~ dunif(0, 10)
+}
+'
+
+# Set up the data
+model_data = list(n = n, y = y, x = x)
+
+# Choose the parameters to watch
+model_parameters =  c("fits")
+
+# Run the model
+model_run = jags(data = model_data,
+                 parameters.to.save = model_parameters,
+                 model.file=textConnection(model_code))
+
+# Simulated results -------------------------------------------------------
+
+# Compute the R-squared value from the fits
+fits = model_run$BUGSoutput$mean$fits
+R_squared = cor(y, fits)^2

practicals/MCMC_game.R

@@ -3,7 +3,7 @@
 # Likelihood is:
 # y[i] ~ normal(alpha, sigma^2)
 # Priors are:
-# alpha ~ normal(0, 2)
+# alpha ~ normal(0, 2^2)
 # sigma ~ uniform(0, 5)
 
 rnorm(1, mean = 0, sd = 2)
@@ -40,10 +40,10 @@ legend('topleft', pch = 1, legend = c('alpha','sigma'), col = c('blue','red'), b
 for(i in 1:iter) {
   alpha_new = rnorm(1, 0, 1)
   sigma_new = runif(1, 0, 5)
-  
+
   post_score_new = sum(dnorm(y, mean = alpha_new, sd = sigma_new, log = TRUE)) + dnorm(alpha, 0, 2, log=TRUE) + dunif(sigma, 0, 5, log=TRUE)
   post_score_old = sum(dnorm(y, mean = alpha, sd = sigma, log = TRUE)) + dnorm(alpha, 0, 2, log=TRUE) + dunif(sigma, 0, 5, log=TRUE)
-  
+
   U = runif(1)
   if(U < exp(post_score_new - post_score_old)) {
     alpha = alpha_new

practicals/tutor_2_R_jags_stan.R

@@ -26,7 +26,7 @@ model {
 }
 '
 model_parameters =  c("intercept", "slope", "residual_sd")
-model_run = jags(data = list(N=nrow(dat),
+model_run = jags(data = list(N = nrow(dat),
                              y = log(dat$earn),
                              x = dat$height_cm),
                  parameters.to.save = model_parameters,

slides/class_1_intro.Rmd

@@ -2,8 +2,9 @@
 title: 'Class 1: An introduction to Bayesian Hierarchical Modelling'
 author: Andrew Parnell \newline \texttt{[email protected]}   \newline \vspace{1cm}
   \newline \includegraphics[width=3cm]{maynooth_uni_logo.jpg}
-  \newline \vspace{2cm}
-  https://andrewcparnell.github.io/bhm_course
+  \newline \vspace{1cm}
+  https://andrewcparnell.github.io/bhm_course 
+  \newline PRESS RECORD 
 output:
   beamer_presentation:
     includes:
@@ -25,7 +26,8 @@ par(mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.01,las=1)
     - where you are from,
     - your previous experience in working with R and regression models,
     - what you are working on,
-    - what you want to get out of the week
+    - what you want to get out of the week,
+    - what you are most looking forward to when we can all go outside again.
 
 - Timetable for the week
 - Pre-requisites
@@ -69,7 +71,7 @@ support your individual studies.
 
 ## Thinking hierarchically: example 1
 
-```{r, echo = FALSE}
+```{r, echo = FALSE, fig.height = 5}
 # Feed conversion rate example
 fcr = read.csv(file = '../data/fcr.csv')
 plot(fcr$food_wt, fcr$weight_gain, 
@@ -79,7 +81,7 @@ plot(fcr$food_wt, fcr$weight_gain,
 
 ## More information:
 
-```{r, echo = FALSE}
+```{r, echo = FALSE, fig.height = 5}
 sex = fcr$sex
 plot(fcr$food_wt[sex=='male'], 
      fcr$weight_gain[sex=='male'], 
@@ -145,7 +147,7 @@ head(swt)
 - The data directory contains a few more data sets which we will play with throughout the week
 - The `data_descriptions.txt` file shows what they contain
 - If you have some spare time it's worth loading them in, exploring relationships, and fitting some simple models
-
+  
 ## Summary
 
 - In hierarchical models we avoid fitting models separately as much as possible

slides/class_2_inference.Rmd

@@ -2,6 +2,7 @@
 title: 'Class 2: Likelihood and Inference'
 author: Andrew Parnell \newline \texttt{[email protected]}   \newline \vspace{1cm}
   \newline \includegraphics[width=3cm]{maynooth_uni_logo.jpg}
+  \newline PRESS RECORD 
 output:
   beamer_presentation:
     includes:
@@ -14,7 +15,7 @@ classoption: "aspectratio=169"
 knitr::opts_chunk$set(dev = 'pdf')
 options(width = 40)
 par(mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.01,las=1)
-pkgs = c('R2jags','rjags', 'lubridate', 'tidyverse','forecast', 'rstan')
+pkgs = c('R2jags', 'rjags', 'lubridate', 'tidyverse', 'rstan')
 lapply(pkgs, library, character.only = TRUE)
 ```
 

slides/class_3_intro_bayes.Rmd

@@ -2,6 +2,7 @@
 title: 'Class 3: An introduction to Bayesian Statistics'
 author: Andrew Parnell \newline \texttt{[email protected]}   \newline \vspace{1cm}
   \newline \includegraphics[width=3cm]{maynooth_uni_logo.jpg}
+  \newline PRESS RECORD 
 output:
   beamer_presentation:
     includes:

slides/class_4_bglms.Rmd

@@ -1,7 +1,8 @@
 ---
-title: 'Class 4: Bayesian Linear and generalised linear models (GLMs)'
+title: 'Class 4: Bayesian linear and generalised linear models (GLMs)'
 author: Andrew Parnell \newline \texttt{[email protected]}   \newline \vspace{1cm}
   \newline \includegraphics[width=3cm]{maynooth_uni_logo.jpg}
+  \newline PRESS RECORD 
 output:
   beamer_presentation:
     includes:

slides/class_6_hglms.Rmd

@@ -2,6 +2,7 @@
 title: 'Class 6: Hierarchical generalised linear models'
 author: Andrew Parnell \newline \texttt{[email protected]}   \newline \vspace{1cm}
   \newline \includegraphics[width=3cm]{maynooth_uni_logo.jpg}
+  \newline PRESS RECORD 
 output:
   beamer_presentation:
     includes:
@@ -142,7 +143,11 @@ library(boot)
 par(mfrow=c(1,3))
 pars = jags_run$BUGSoutput$mean
 for(i in c(2,3,1)) {
-  plot(swt$forest, y/N, main = paste('Altitude type:',levels(swt$alt)[i]), xlab = '% forest cover', ylab = 'Estimated proporton')
+  curr_rows = which(swt$alt == levels(swt$alt)[i])
+  plot(swt$forest[curr_rows], y[curr_rows]/N[curr_rows], 
+       main = paste('Altitude type:',levels(swt$alt)[i]),
+       xlab = '% forest cover', 
+       ylab = 'Estimated proporton')
   points(swt$forest, inv.logit(pars$alpha[i] + pars$beta[i] * swt$forest), col = i+1)
 }
 par(mfrow=c(1,1))

slides/class_7_multihms.Rmd

@@ -2,6 +2,7 @@
 title: 'Class 7: Multi-layer hierarchical models'
 author: Andrew Parnell \newline \texttt{[email protected]}   \newline \vspace{1cm}
   \newline \includegraphics[width=3cm]{maynooth_uni_logo.jpg}
+  \newline PRESS RECORD 
 output:
   beamer_presentation:
     includes:

slides/class_8_pp_zi_multinom.Rmd

@@ -2,6 +2,7 @@
 title: 'Class 8: Partial pooling and zero-inflation'
 author: Andrew Parnell \newline \texttt{[email protected]}   \newline \vspace{1cm}
   \newline \includegraphics[width=3cm]{maynooth_uni_logo.jpg}
+  \newline PRESS RECORD 
 output:
   beamer_presentation:
     includes:

slides/class_9_varsel_JAGS.Rmd

@@ -2,6 +2,7 @@
 title: 'Class 9: Variable selection and non-parametrics'
 author: Andrew Parnell \newline \texttt{[email protected]}   \newline \vspace{1cm}
   \newline \includegraphics[width=3cm]{maynooth_uni_logo.jpg}
+  \newline PRESS RECORD 
 output:
   beamer_presentation:
     includes: