Nothing
#' Prior Predictive Replication p-value Calculation
#'
#' Assessing the prior predictive distribution and calculating the replication
#' p-value based on it.
#'
#' @param beta A 2-D vector, containing the estimates in the original study and the
#' replication study.
#' @param se A 2-D vector, containing the standard errors of the estimates in the original
#' study and the replication study.
#' @param r_vec A vector, defining the prior reproducible model. Each r value
#' corresponds to a probability of sign consistency.
#' @param test A string, determining which test statistics to utilize. If not specified,
#' the default two-sided one will be used.
#' @param report_PI A boolean, denoting whether the 95% predictive interval for
#' the estimates be reported or not. This option is only valid for two-sided
#' test statistics. The default is FALSE.
#'
#' @return
#' A list with the following components:
#' \item{grid}{ The detailed grid values for the hyperparameters.}
#' \item{test_statistics}{ The test statistics used in calculating the replication p-value.}
#' \item{pvalue}{ The resulting prior predictive replicaiton p-value.}
#' \item{predictive_interval}{The 95% predictive interval if required.}
#'
#' @importFrom stats dnorm integrate pnorm qchisq rnorm uniroot
#'
#' @export
#'
#' @examples
#' data("RPP_filtered")
#' attach(RPP_filtered)
#' rpp_pval<-sapply(1:nrow(RPP_filtered),function(x)
#' prior_prp(beta=c(beta_orig[x], beta_rep[x]),se=c(se_orig[x], se_rep[x]))$pvalue)
#'
prior_prp <- function(beta,se,r_vec = c(0, 8e-4, 6e-3, 0.024),test="two_sided",report_PI=FALSE){
reslist <- list()
beta_o <- beta[1]
beta_r <- beta[2]
se2_o <- se[1]^2
se2_r <- se[2]^2
## grid set-up
#chis1<-qchisq(c(0.25,0.5,0.75),df=1)
eta2_vec = (se2_o+beta_o^2)/qchisq(c(0.25, 0.5, 0.75), df=1)
#pv = c(1.0, 0.99, 0.975, 0.95)
rv = r_vec
make_grid <-function(eta2){
grid = sapply(rv, function(x) c(eta2*(1-x), eta2*x))
return(t(grid))
}
grid = c()
for (i in 1:length(eta2_vec)){
grid = rbind(grid, make_grid(eta2_vec[i]))
}
reslist[["grid"]] = grid
omg2 = grid[,1]
phi2 = grid[,2]
# prior updates
mean <- sapply(1:length(omg2),function(x)
beta_o/(se2_o/omg2[x]+phi2[x]/omg2[x]+1))
var <- sapply(1:length(omg2),function(x)
(se2_o+phi2[x])*omg2[x]/(se2_o+phi2[x]+omg2[x])+phi2[x]+se2_r)
pval <- sapply(1:length(mean),function(x)
pnorm(beta_r, mean=mean[x],sd=sqrt(var[x])))
#weights for each model
wts = dnorm(beta_o, mean=0, sd=sqrt(omg2+se2_o+phi2))
wts = wts/sum(wts)
pval_wt <- wts%*%pval
res = NA
reslist["test_statistics"] = test
if (test=="pub_bias"){
mean_scale <- mean/beta_o
sd_scale <- sqrt(var)/abs(beta_o)
pval_new <- sapply(1:length(mean),function(x)
pnorm(beta_r/beta_o, mean=mean_scale[x],sd=sd_scale[x]))
res = wts%*%pval_new
reslist[["pvalue"]] = res
}
if (test == "two_sided"){
res <- 2*min(pval_wt,1-pval_wt)
reslist[["pvalue"]]=res
if (report_PI){
mixture_CDF_right <- Vectorize(function(t)
wts%*%pnorm(t,mean=mean,sd=sqrt(var))-0.975)
root_right <- uniroot(mixture_CDF_right,c(mean(mean),10))$root
mixture_CDF_left <- Vectorize(function(t)
wts%*%pnorm(t,mean=mean,sd=sqrt(var))-0.025)
root_left <- uniroot(mixture_CDF_left,c(-10,mean(mean)))$root
reslist[["predictive_interval"]] <- c(PI_left_95p=root_left,PI_right_95p=root_right)
}
}
return(reslist)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.