Spatial-Temporal Priors in the STAP model

library(rstap)
library(dplyr)
library(tidyr)
library(ggplot2)
library(purrr)
library(plot3D)

Starting with a simple model

\[ E[Y] = \alpha + Z_1\delta_{sex} + X(\theta)_{FF}\beta_{FF} \] We’ll calculate the likelihood for different values of \(\theta\) and \(\beta_FF\) assuming that \(\alpha and \delta\) are given

data("homog_subject_data")
data("homog_distance_data")
homog_subject_data
homog_distance_data
alpha <- 22
delta <- -.8
sigma <- 2.3
thetas = seq(from = 0.01, to = 10, by =0.1)
betas = seq(from = 0.01, to = 10, by = 0.1)
Z <- model.matrix(y ~ sex, data = homog_subject_data)
Xs <- homog_distance_data %>% crossing(thetas)  %>% group_by(subj_id,thetas) %>% 
    summarise(exposure = sum(pracma::erfc(Distance/thetas))) %>% arrange(thetas) %>% 
    split(.$thetas)
Xs <- map(Xs,function(x) x %>% mutate(X_tilde = (exposure - mean(x$exposure))/sd(x$exposure)))
lls <- map(betas,function(a) map_dbl(Xs,function(x)  sum(dnorm(x = homog_subject_data$y,mean = Z %*% c(alpha,delta) +   x$X_tilde * a,sd = sigma,log = T))))
llmat <- Reduce(f = rbind,x = lls)
rownames(llmat) <- as.character(betas)
colnames(llmat) <- as.character(thetas)
par(mfrow=c(1,2))
persp3D(x = betas, y = thetas, z = llmat,phi = .6,theta = 0.5,xlab="beta",ylab="theta",zlab="log-lik",ticktype ='detailed')
persp3D(x = betas, y = thetas, z = llmat, phi = 3, theta = 90, xlab="beta",ylab="theta",zlab="log-lik",ticktype = "detailed")
data_frame(beta = betas[1],
           thetas = thetas,
           ll = llmat[1,]) %>% 
    rbind(data_frame(beta = betas[15],
                     thetas = thetas,
                     ll = llmat[15,])) %>% 
    rbind(data_frame(beta = betas[50],
                     thetas = thetas,
                     ll = llmat[50,])) %>% 
    rbind(data_frame(beta = betas[75],
                     thetas = thetas,
                     ll = llmat[75,])) %>% 
    rbind(data_frame(beta = betas[100],
                     thetas = thetas,
                     ll = llmat[100,])) %>% 
    mutate(beta = factor(beta)) %>% 
    ggplot(aes(x=thetas, y= ll, color = beta )) + geom_line() + theme_bw() + 
    geom_vline(aes(xintercept = .5),linetype = 2, color= 'black') + 
    ggtitle("Log - likelihood for possible values of theta, beta")
data_frame(beta = betas,
           theta = thetas[1],
           ll = llmat[,1]) %>% 
    rbind(data_frame(beta = betas,
                     theta = thetas[25],
                     ll = llmat[,25])) %>% 
    rbind(data_frame(beta = betas,
                     theta = thetas[50],
                     ll = llmat[,50])) %>% 
    rbind(data_frame(beta = betas,
                     theta = thetas[100],
                     ll = llmat[,100])) %>% 
    mutate(theta = factor(theta)) %>% 
    ggplot(aes(x=beta, y= ll, color = theta )) + geom_line() + theme_bw() + 
    geom_vline(aes(xintercept = 1.2),linetype = 2, color= 'black') + 
    ggtitle("Log - likelihood for possible values of beta, theta")

\[ p(\theta,\beta|Y,\alpha,\delta,\sigma) \propto p(y|...)p(\beta)p(\theta)\\ \beta \sim N(0,3)\\ \log(\theta) \sim N(1,1) \]

post <- map(1:length(betas),function(a) map_dbl(Xs,function(x)  sum(dnorm(x = homog_subject_data$y,mean = Z %*% c(alpha,delta) +   x$X_tilde * betas[a],
                                                                             sd = sigma,log = F)) * 
                                  dnorm(x = betas[a],mean = 0,sd = 3,log = F) * 1)*
                                  dlnorm(x = thetas[a],meanlog = 1, sdlog = 1, log =F))
postmat <- Reduce(f = rbind,x = post)
par(mfrow=c(1,2))
persp3D(x = betas, y = thetas, z = postmat, phi = .6,theta = 0.5,xlab="beta",ylab="theta", zlab="density", ticktype = "detailed", main= "Posterior-unnormalized")
persp3D(x = betas, y = thetas, z = postmat, phi = 3, theta = 90, xlab="beta",ylab="theta", zlab="density", ticktype = "detailed")

Looking at the \(\theta\) cross section, specifically

data_frame(beta = betas[1],
           thetas = thetas,
           post_dens = postmat[1,]) %>% 
    rbind(data_frame(beta = betas[50],
                     thetas = thetas,
                     post_dens = postmat[50,])) %>% 
    rbind(data_frame(beta = betas[100],
                     thetas = thetas,
                     post_dens = postmat[100,])) %>% 
    rbind(data_frame(beta = betas[120],
                     thetas = thetas,
                     post_dens = postmat[120,])) %>% 
    rbind(data_frame(beta = betas[150],
                     thetas = thetas,
                     post_dens = postmat[150,])) %>% 
    mutate(beta = factor(beta)) %>% 
    ggplot(aes(x=thetas, y= post_dens, color = beta )) + geom_line() + theme_bw() + 
    geom_vline(aes(xintercept = .5),linetype = 2, color= 'black') + 
    ggtitle("Posterior Density (unnormalized) for possible values of theta, beta")