Skip to content

Commit

Permalink
start wrapper
Browse files Browse the repository at this point in the history
  • Loading branch information
@s-huebler
s-huebler committed Oct 17, 2024
1 parent 5b9b473 commit 6b19208
Show file tree
Hide file tree
Showing 5 changed files with 124 additions and 98 deletions.
73 changes: 49 additions & 24 deletions Functions/gibbs.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,37 +6,62 @@
#
# Purpose: Single iteration of a gibbs sampling
# Advanced R Goals:
# 1) Remove for loop for the simulations (parallelized later)
# 1) Remove for loop for the simulations
# 2) Modularization
# Notes:
# - Use set_hyperparameters function to define hypers
# - Source mh function
#############################################################'

#Gibbs implementation - update mu, tau2, gammas, thetas
gibbs <- function(mu,
tau2,
gammas,
thetas,
k,
hypers){
gibbs <- function(Y, #Data
mu, #Current mu estimate
tau2, #Current tau2 estimate
gammas, #Vector size k
thetas, #Vector size k
k, #Number of observations in meta-analysis
hypers #Hyperparameters
){

#Update mu first
mu_new <- post_mu(k=k,
sigma = 1/tau2,
theta = thetas,
b_mu = hypers$b_mu)


#Update tau2 using the mu_new
tau2_new <- post_tau(k = k,
a_tau = hypers$a_tau2,
b_tau = hypers$b_tau2,
mu = mu_new,
theta = thetas)



#Update mu first
mu_new <- post_mu(k=k,
sigma = 1/tau2,
theta = thetas,
b_mu = hypers$b_mu)


#Update tau2 using the mu_new
tau2_new <- post_tau(k = k,
a_tau = hypers$a_tau2,
b_tau = hypers$b_tau2,
mu = mu_new,
theta = thetas)

mh_iter <- mh(Y,
mu_new,
tau2_new,
gammas,
thetas,
k,
hypers

)

ret <- list("mu" = mu_new,
"tau2" = tau2_new,
"thetas" = mh_iter$thetas,
"gammas" = mh_iter$gammas)

ret

}


}
gibbs(dat,
mu_test,
tau2_test,
gammas_test,
thetas_test,
7,
hypers
)
15 changes: 15 additions & 0 deletions Functions/gibbs_mh_run.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,15 @@
############################################################'
# gibbs_mh_run.R
#
# Author: Sophie Huebler
# Creation Date: 15Oct2024
#
# Purpose: Wrapper for full algorithm
#
# Advanced R Goals:
# 1) Use data table
#
#
# Notes: - Use set_hyperparameters function to define hypers
#############################################################'

77 changes: 32 additions & 45 deletions Functions/mh.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,6 @@
# 1) Remove lapply statements replace with single loop
# (maybe switch to cpp?)
# 2) Modularization
# 3) Pre allocating storage
#
#
# Notes: - Use set_hyperparameters function to define hypers
Expand All @@ -32,17 +31,11 @@ mh <- function(Y, #Data
hypers #Hyperparameters
){

#Pre allocate space for vectors from for loop
##Current posteriors
current_post_thetas <- numeric(k)
current_post_gammas <- numeric(k)

##Candidate posteriors
candidate_post_thetas <- numeric(k)
candidate_post_gammas <- numeric(k)


#Pre allocate memory
thetas_new <- thetas
gammas_new <- gammas


#Define candidates
candidates_thetas <- rnorm(k,
mean(as.numeric(thetas)),
Expand All @@ -52,72 +45,66 @@ mh <- function(Y, #Data
mean(as.numeric(gammas)),
hypers$mh_scale)

#Define critical value
crit <- log(runif(1))

#Run posterior function on current and proposed data
# for thetas and gammas in a single loop
# for thetas and gammas and compare to crit value in a single loop

for(i in 1:k){

#Thetas
current_post_thetas[i] <- post_theta(tau2 = tau2,
temp_theta_current <- post_theta(tau2 = tau2,
mu = mu,
gamma = gammas[i],
theta = thetas[i],
rt = Y[i, "rit"],
rc = Y[i, "ric"]) |>
log()

candidate_post_thetas[i] <- post_theta(tau2 = tau2,
temp_theta_candidate <- post_theta(tau2 = tau2,
mu = mu,
gamma = candidates_gammas[i],
theta = candidates_thetas[i],
rt = Y[i, "rit"],
rc = Y[i, "ric"]) |>
log()

if(crit < (temp_theta_candidate - temp_theta_current)){
thetas_new[i]<- candidates_thetas[i]
}

#Gammas
current_post_gammas[i] <- post_gamma(tau2 = tau2,
temp_gamma_current <- post_gamma(b_gamma = hypers$b_gamma,
tau2 = tau2,
mu = mu,
gamma = gammas[i],
theta = thetas[i],
rt = Y[i, "rit"],
rc = Y[i, "ric"]) |>
log()

candidate_post_gammas[i] <- post_gamma(tau2 = tau2,
temp_gamma_candidate <- post_gamma(b_gamma = hypers$b_gamma,
tau2 = tau2,
mu = mu,
gamma = candidates_gammas[i],
theta = candidates_thetas[i],
rt = Y[i, "rit"],
rc = Y[i, "ric"]) |>
log()

}

#Divide posteriors = subtract log posteriors
logr_theta = proposed_post_theta - current_post_theta
logr_gamma = proposed_post_gamma - current_post_gamma

#Update based R function
#Using uniformdensity
crit <- log(runif(1))
criteria_theta <- crit < logr_theta
criteria_gamma <- crit < logr_gamma

#Vectorization

#Inner loop (fix later for efficiency)
#Update estimates if ratio of posterior of current
#and proposed estimates is less than the log function
for(i in 1:k){
if(criteria_theta[i]){
thetas[j,i]<-candidates_thetas[i]
}else{
thetas[j,i] <- thetas[(j-1),i]
}
if(criteria_gamma[i]){
gammas[j,i]<-candidates_gammas[i]
}else{
gammas[j,i] <- gammas[(j-1),i]
if(crit < (temp_gamma_candidate - temp_gamma_current)){
gammas_new[i]<- candidates_gammas[i]
}
}
}

rm(temp_gamma_candidate,
temp_gamma_current,
temp_theta_candidate,
temp_theta_current)
} #End for loop

ret <- list("thetas" = thetas_new,
"gammas" = gammas_new)

ret
}
54 changes: 27 additions & 27 deletions Original_Code.R
View file
Original file line number Diff line number Diff line change
@@ -1,36 +1,36 @@
#functions based on conditional posteriors

#closed form for mu and tau
post_mu <- function(k, sigma, theta, b_mu){
rnorm(1,
mean = sigma*sum(theta)/(k*sigma + 100),
sd = 1/sqrt(k*sigma + 100))
}

post_tau <- function(k, a_tau, b_tau, mu, theta){
1/rgamma(1,
shape = k/2 + a_tau,
rate = 0.5*sum((theta-mu)^2) + b_tau)
}

#no closed form for gammas and thetas
post_theta <- function(tau2, mu, gamma, theta, rt, rc){
exp(theta/2*(rt-rc))/
((1+exp(gamma - theta/2))^(rc)*(1+exp(gamma + theta/2))^(rt))*
1/sqrt(2*pi*tau2)*
exp(-1/(2*tau2)*(theta-mu)^2)
}

post_gamma <- function(b_gamma, mu, tau2, theta, gamma, rt, rc){
exp(gamma*(rt-rc))/
((1+exp(gamma - theta/2))^(rc)*(1+exp(gamma + theta/2))^(rt))*
1/sqrt(2*pi*b_gamma)*
exp(-1/(2*tau2)*(gamma)^2)
}
# post_mu <- function(k, sigma, theta, b_mu){
# rnorm(1,
# mean = sigma*sum(theta)/(k*sigma + 100),
# sd = 1/sqrt(k*sigma + 100))
# }
#
# post_tau <- function(k, a_tau, b_tau, mu, theta){
# 1/rgamma(1,
# shape = k/2 + a_tau,
# rate = 0.5*sum((theta-mu)^2) + b_tau)
# }
#
# #no closed form for gammas and thetas
# post_theta <- function(tau2, mu, gamma, theta, rt, rc){
# exp(theta/2*(rt-rc))/
# ((1+exp(gamma - theta/2))^(rc)*(1+exp(gamma + theta/2))^(rt))*
# 1/sqrt(2*pi*tau2)*
# exp(-1/(2*tau2)*(theta-mu)^2)
# }
#
# post_gamma <- function(b_gamma, mu, tau2, theta, gamma, rt, rc){
# exp(gamma*(rt-rc))/
# ((1+exp(gamma - theta/2))^(rc)*(1+exp(gamma + theta/2))^(rt))*
# 1/sqrt(2*pi*b_gamma)*
# exp(-1/(2*tau2)*(gamma)^2)
# }



mh_gibbs <- function(#Y a data frame with rit, nit, ric, nic, theta_i, gamma_i
mh_gibbs_original <- function(#Y a data frame with rit, nit, ric, nic, theta_i, gamma_i
Y,

#hyperparams
Expand Down
3 changes: 1 addition & 2 deletions data/read_data.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,5 +42,4 @@ dat <- metafor::escalc(measure = "OR",
dplyr::mutate(sd_i = sqrt(v_i))|>
dplyr::mutate(lower_ci = theta_i - 1.96*sd_i,
upper_ci = theta_i + 1.96*sd_i)|>
dplyr::mutate(gamma_i = 0.5*(qlogis(rit/nit)+qlogis(ric/nic)))|>
dplyr::select(trial:v_i)
dplyr::mutate(gamma_i = 0.5*(qlogis(rit/nit)+qlogis(ric/nic)))

0 comments on commit 6b19208

Please sign in to comment.