v6_li_mat_func <- function(pars, data, grid, type=c("null", "causal", "shared")){
type <- match.arg(type)
rho <- pars[1]
if(type %in% c("causal", "shared")){
stopifnot(length(pars) > 1)
b <- pars[2]
}else{
b <- 0
}
ll_mat <- ll_v6_mat(rho, b,
grid[,1],
data$beta_hat_1,
data$beta_hat_2,
data$seb1,
data$seb2)
cols <- c(0, rep(1:3, each=nrow(grid)))
if(type=="null"){
ll_mat <- ll_mat[, cols!=2]
}else if(type=="causal"){
ll_mat <- ll_mat[, cols!=1]
}
return(ll_mat)
}
v6_prior_func <- function(pars, type=c("null", "causal", "shared"),
z_prior_func=NULL,
b_prior_func = NULL){
type <- match.arg(type)
arctanh <- function(rho){
0.5*log((1+rho)/(1-rho))
}
if(is.null(z_prior_func)){
z_prior_func = function(z){dnorm(z, 0, 0.5, log=TRUE)}
}
if(is.null(b_prior_func)){
b_prior_func <- function(b){
log((b/0.4)^2) + dnorm(b, 0, 0.4, log=TRUE)
}
}
rho <- pars[1]
z <- arctanh(rho)
ll <- z_prior_func(z)
if(type %in% c("causal", "shared")){
stopifnot(length(pars) > 1)
b <- pars[2]
ll <- ll + b_prior_func(b)
}
return(ll)
}
v6_q_gamma_from_pi <- function(pi, type=c("null", "causal", "shared")){
n <- length(pi)
if(type=="null"){
L <- (n-1)/2
stopifnot(round(L)==L)
gamma <- c(pi[1], sum(pi[2:(L+1)]),
0,sum(pi[(L + 2):n]))
q <- gamma[3]/(sum(gamma[2:3]))
pi1 <- pi[2:(L+1)]/gamma[2]
pi2 <- rep(NA, L)
pi3 <- pi[(L + 2):n]/gamma[4]
return(list("gamma"=gamma, "q"=q,
"pi1"=pi1, "pi2"=pi2, "pi3"=pi3))
}
if(type=="causal"){
L <- (n-1)/2
stopifnot(round(L)==L)
gamma <- c(pi[1], 0, sum(pi[2:(L+1)]),
sum(pi[(L + 2):n]))
q <- gamma[3]/(sum(gamma[2:3]))
pi1 <- rep(NA, L)
pi2 <- pi[2:(L+1)]/gamma[3]
pi3 <- pi[(L + 2):n]/gamma[4]
return(list("gamma"=gamma, "q"=q,
"pi1"=pi1, "pi2"=pi2, "pi3"=pi3))
}
if(type=="shared"){
L <- (n-1)/3
stopifnot(round(L)==L)
gamma <- c(pi[1], sum(pi[2:(L+1)]),
sum(pi[(L+2):(2*L + 1)]),
sum(pi[(2*L + 2):n]))
q <- gamma[3]/(sum(gamma[2:3]))
pi1 <- pi[2:(L+1)]/gamma[2]
pi2 <- pi[(L+2):(2*L + 1)]/gamma[3]
pi3 <- pi[(2*L + 2):n]/gamma[4]
return(list("gamma"=gamma, "q"=q,
"pi1"=pi1, "pi2"=pi2, "pi3"=pi3))
}
}
v6_li_func_share <- function(pars, data, grid, gamma_red,
z_prior_func, b_prior_func){
#pars: z, b, logitq
z <- pars[1]
rho <- tanh(z)
b <- pars[2]
logitq <- pars[3]
q <- expit(pars[3])
#grid should have first column as non-zero sigma and next three columns as pi1, pi2, pi3
gamma <- c(gamma_red[1], (1-q)*gamma_red[2], q*gamma_red[2], gamma_red[3])
pi <- c(gamma[1], gamma[2]*grid[,2], gamma[3]*grid[,3], gamma[4]*grid[,4])
K <- length(pi)
ll1 <- ll_v6(rho, b, gamma,
grid[,2], grid[,3], grid[,4], grid[,1],
data$beta_hat_1, data$beta_hat_2,
data$seb1, data$seb2, data$wts)
prior <- z_prior_func(z) + b_prior_func(b) + ddirichlet1(pi, rep(c(10, 1), c(1, K-1)))
return(ll1 + prior)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.