Rcode.r

library(ggplot2)
library(tictoc)

######################################################
################ The usual approach ##################
######################################################

load("~/OXTsum.RData")
OXTsum$Species <- as.factor(OXTsum$Species)
OXTsum$Temperature <- as.factor(OXTsum$Temperature)
OXTsum$TotalArea <- as.numeric(OXTsum$TotalArea)

Boxplots <- ggplot(OXTsum, aes(x=Temperature, y=TotalArea, colour=Species))
Boxplots + geom_boxplot()

# Difference in slopes? Hard with vole scatter
load("~/OXTsum.RData")
OXTsum$Species <- as.factor(OXTsum$Species)
OXTsum$Sex <- as.factor(OXTsum$Sex)
OXTsum$Temperature <- as.numeric(OXTsum$Temperature)
OXTsum$TotalArea <- as.numeric(OXTsum$TotalArea)

linreg <- lm(TotalArea ~ Temperature*Species, data = OXTsum)
summary(linreg)

# Just Mouse slope:
justmice <- OXTsum[OXTsum$Species == 'Mouse',]
linreg0 <- lm(TotalArea ~ 1, data = justmice)
linreg1 <- lm(TotalArea ~ Temperature, data = justmice)
linreg2 <- lm(TotalArea ~ Temperature+Sex, data = justmice)
summary(linreg1)
anova(linreg0, linreg1, linreg2) # p=0.04

# Just Vole slope:
justvoles <- OXTsum[OXTsum$Species == 'Vole',]
linreg0 <- lm(TotalArea ~ 1, data = justvoles)
linreg1 <- lm(TotalArea ~ Temperature, data = justvoles)
linreg2 <- lm(TotalArea ~ Temperature+Sex, data = justvoles)
summary(linreg1)
anova(linreg0, linreg1, linreg2) # p=0.695


######################################################
################# Bayesian approach ##################
######################################################


loglikelihood <- function(iV, iM, sV, sM, slopeV, slopeM){
  ll <- 0
  for (idx in 1:nrow(OXTsum)){
    x           <- OXTsum$TotalArea[idx]
    species     <- OXTsum$Species[idx]
    temperature <- OXTsum$Temperature[idx]
    if (species == 'Mouse'){
      ll = ll  - ( (x-(iM + slopeM*temperature))^2 / (2*sM*sM) ) - 0.5*log(2*pi*sM*sM)
    } else {
      ll = ll - ( (x-(iV + slopeV*temperature))^2 / (2*sV*sV) ) - 0.5*log(2*pi*sV*sV)
    }
  }
  return(ll)
}

negll <- function(x){
  iV  <- x[1]
  iM  <- x[2]
  sV  <- x[3]
  sM <- x[4]
  slopeV <- x[5]
  slopeM <- x[6]
  return( -1* loglikelihood( iV, iM, sV, sM, slopeV, slopeM) )
}

iM <- mean(OXTsum[OXTsum$Species == 'Mouse',]$TotalArea)
iV <- mean(OXTsum[OXTsum$Species == 'Vole',]$TotalArea)
sM <- sd(OXTsum[OXTsum$Species == 'Mouse',]$TotalArea)
sV <- sd(OXTsum[OXTsum$Species == 'Vole',]$TotalArea)
slopeM <- 0
slopeV <- 0

MAP <- optim(c(iV, iM, sV, sM, slopeV, slopeM), negll) # "Maximum a posteriori"
# 131645.2285  35317.9663  71994.4467   9742.3578    764.0969   -683.8436

# Mouse slope different from zero? -- Laplace approximation. Fast
H <- numDeriv::hessian(func = negll, x = MAP$par)
H <- (H + t(H)) / 2 #  make symmetric
cov_mat <- (solve(H) + t(solve(H)))/2
mean_mouse <- MAP$par[6]
sd_mouse  <- sqrt(cov_mat[6, 6])
mean_mouse/sd_mouse

alpha <- 0.05
ci_lower <- mean_mouse + qnorm(alpha/2) * sd_mouse
ci_upper <- mean_mouse + qnorm(1 - alpha/2) * sd_mouse
pnorm(mean_mouse/sd_mouse) # values below z-score


# Mouse slope different from zero? -- Monte Carlo Sampling to Marginalize

marginalize_mouse <- function(negll, MouseSlopeValue, nsamp = 10000) {
  idx_marg <- 1:5   
  mu_marg <- MAP$par[idx_marg]
  cov_marg <- cov_mat[idx_marg, idx_marg]
  samples_marg <- MASS::mvrnorm(nsamp, mu_marg, cov_marg)
  MouseSlope_rep <- rep(MouseSlopeValue, each = nsamp)
  param_rep <- samples_marg[rep(1:nsamp, times = length(MouseSlopeValue)),,drop = FALSE]
  X_all <- cbind(param_rep, MouseSlope = MouseSlope_rep)
  loglik_all <- apply(X_all, 1, function(x) -negll(x))   # Evaluate log-likelihoods for all (samples, MouseSlope) combinations
  loglik_mat <- matrix(loglik_all, nrow = nsamp, ncol = length(MouseSlopeValue))
  max_loglik <- apply(loglik_mat, 2, max)
  marg_ll <- log(colMeans(exp(loglik_mat - matrix(max_loglik, nrow = nsamp, ncol = length(MouseSlopeValue), byrow = TRUE)))) + max_loglik
  names(marg_ll) <- MouseSlopeValue
  marg_ll
}

tic()
mouseSlopeValues <- seq(-2000, 500, length.out = 100)
result_1D <- marginalize_mouse(negll, mouseSlopeValues) # Will take 20-30 min
toc()

pp <- exp(result_1D - max(result_1D))
pp <- pp / (sum(pp)*mean(diff(mouseSlopeValues)))  # Normalize for proper posterior
plot(mouseSlopeValues, pp, type = "l", xlab = "Mouse Slope", ylab = "Posterior Likelihood")
sum(pp[mouseSlopeValues>=0]) / sum(pp)


# Mouse + Vole  - 2D Sampling
marginalize_mouse_vole <- function(negll, MouseSlopeValues, VoleSlopeValues, nsamp = 10000) {
  idx_marg <- 1:4
  mu_marg <- MAP$par[idx_marg]
  cov_marg <- cov_mat[idx_marg, idx_marg]
  samples_marg <- MASS::mvrnorm(nsamp, mu_marg, cov_marg)
  grid <- expand.grid(MouseSlope = MouseSlopeValues, VoleSlope = VoleSlopeValues)
  n_grid <- nrow(grid)
  param_rep <- samples_marg[rep(1:nsamp, times = n_grid), , drop = FALSE]
  grid_rep <- grid[rep(1:n_grid, each = nsamp), , drop = FALSE]
  X_all <- cbind(param_rep, VoleSlope = grid_rep$VoleSlope, MouseSlope = grid_rep$MouseSlope)
  loglik_all <- apply(X_all, 1, function(x) -negll(x))   # Evaluate log-likelihood for all points
  loglik_mat <- matrix(loglik_all, nrow = nsamp, ncol = n_grid) # Reshape results into nsamp × n_grid matrix
  max_loglik <- apply(loglik_mat, 2, max)
  marg_ll_vec <- log(colMeans(exp(loglik_mat - matrix(max_loglik, nrow = nsamp, ncol = n_grid, byrow = TRUE)))) + max_loglik # Compute marginal log-likelihood (log-sum-exp trick)
  marg_ll <- matrix(marg_ll_vec, nrow = length(MouseSlopeValues), ncol = length(VoleSlopeValues), dimnames = list(MouseSlope = MouseSlopeValues, VoleSlope = VoleSlopeValues))
  marg_ll
}


MouseSlopeValues <- seq(-2000, 500, length.out = 100)
VoleSlopeValues <- seq(-5000, 6000, length.out = 100)

tic()
result_2D <- marginalize_mouse_vole(negll, MouseSlopeValues, VoleSlopeValues)
toc()

image(MouseSlopeValues, VoleSlopeValues, exp(result_2D - max(result_2D)), xlab = "Mouse Slope", ylab = "Vole Slope", main = "2D Likelihood")