From cf8a996cf8349073d5bcd8df955a3d9008c921b9 Mon Sep 17 00:00:00 2001 From: Sophie Huebler Date: Tue, 15 Oct 2024 11:46:36 -0600 Subject: [PATCH] Initial commit from rstudio --- .gitignore | 3 +- Original_Code.R | 200 ++++++++++++++++++++++++++++++++++++++++++ PHS7045_Midterm.Rproj | 13 +++ 3 files changed, 215 insertions(+), 1 deletion(-) create mode 100644 Original_Code.R create mode 100644 PHS7045_Midterm.Rproj diff --git a/.gitignore b/.gitignore index b7597b3..ac6ed54 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ *_cache -*_files \ No newline at end of file +*_files +.Rproj.user diff --git a/Original_Code.R b/Original_Code.R new file mode 100644 index 0000000..e3af948 --- /dev/null +++ b/Original_Code.R @@ -0,0 +1,200 @@ +#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) +} + + + +mh_gibbs <- function(#Y a data frame with rit, nit, ric, nic, theta_i, gamma_i + Y, + + #hyperparams + #mean and var for mu + a_mu = 0, + b_mu = 100, + #mean and var for gamma + a_gamma = 0, + b_gamma = 100, + #shape and scale for 1/tau2 + a_tau2 = 0.001, + b_tau2 = 0.001, + + #Number of simulations + sim = 1e5, + + #Burn + burn = 2000){ + + # Set up + #number of studies + k <- nrow(Y) + + #initial values based on data and hyperparameters + #thetas initialized by odds ratio in data + thetas <- data.frame(matrix(as.vector(Y[,"theta_i"]), + nrow = 1))%>% + mutate_all(as.numeric) + + #gammas initialized by average logit in data + gammas <- data.frame(matrix(as.vector(Y[,"gamma_i"]), + nrow = 1))%>% + mutate_all(as.numeric) + + #mu initialized by average theta + mu <- mean(thetas[1,]%>%as.numeric()) + + #tau initialized by variance of thetas + tau2 <- sd(thetas)^2 + + + #Outer loop - number of simulations + for(j in 2:sim){ + #Gibbs implementation - start with mu,tau + + mu_new <- post_mu(k=k, + sigma = 1/tau2[(j-1)], + theta = thetas[(j-1),], + b_mu = b_mu) + mu <- c(mu, mu_new) + + tau2_new <- post_tau(k = k, + a_tau = a_tau2, + b_tau = b_tau2, + mu = mu[j], + theta = thetas[(j-1),]) + + tau2 <- c(tau2, tau2_new) + + + #Metropolis-Hastings procedure for + #componentwise updating of theta, gamma + + #Candidate density: Normal centered on + #current value and variance low (suggested + #by Dan in class for scaling 0.23) + + candidates_thetas <- rnorm(k, mean(thetas[(j-1),]%>% + as.numeric()), 0.23) + candidates_gammas <- rnorm(k, mean(gammas[(j-1),]%>% + as.numeric()), 0.23) + + ##Run posterior function on current and proposed data + + #Theta + #Find the log posterior for the currrent estimates + current_post_theta <- lapply(1:k, function(i) + post_theta( + tau2 = tau2[j], + mu = mu[j], + gamma = gammas[(j-1),i], + theta = thetas[(j-1), i], + rt = Y[i, "rit"], + rc = Y[i, "ric"]) + ) %>% + unlist() %>% + log() + + #Find the log posterior for the proposed estimates + proposed_post_theta <- lapply(1:k, function(i) + post_theta( + tau2 = tau2[j], + mu = mu[j], + gamma = candidates_gammas, + theta = candidates_thetas, + rt = Y[i, "rit"], + rc = Y[i, "ric"])) %>% + unlist() %>% + log() + + #Gamma + #Find the log posterior for the current estimates + current_post_gamma<- lapply(1:k, function(i) + post_gamma( + b_gamma = b_gamma, + mu = mu[j], + tau2 = tau2[j], + gamma = gammas[(j-1),i], + theta = thetas[(j-1), i], + rt = Y[i, "rit"], + rc = Y[i, "ric"]) + ) %>% + unlist() %>% + log() + + #Find the log posterior for the proposed estimates + proposed_post_gamma <- lapply(1:k, function(i) + post_gamma( + b_gamma = b_gamma, + mu = mu[j], + tau2 = tau2[j], + gamma = candidates_gammas, + theta = candidates_thetas, + rt = Y[i, "rit"], + rc = Y[i, "ric"])) %>% + unlist() %>% + 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 + criteria_theta <- log(runif(1)) < logr_theta + criteria_gamma <- log(runif(1)) < logr_gamma + + #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] + } + }#end inner loop + + }#end outer loop + + #print(thetas[(burn+1:nrow(thetas)),]) + + #Burn + ret <-list("mu" = mu[(burn+1):length(mu)], + "tau2" = tau2[(burn+1):length(tau2)], + "thetas" = thetas[(burn+1:nrow(thetas)),]%>% + filter(complete.cases(.)), + "gammas" = gammas[(burn+1:nrow(gammas)),]%>% + filter(complete.cases(.))) + ret +} diff --git a/PHS7045_Midterm.Rproj b/PHS7045_Midterm.Rproj new file mode 100644 index 0000000..8e3c2eb --- /dev/null +++ b/PHS7045_Midterm.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX