# Chapter 7 Supplementary Material

## 7.1 JAGS Models

cond.likelihood.model = function() {
for (j in 1:J) {
CaseCon[j] ~ dbern(theta[j])
logit(theta[j]) <-  mu.p53  }

C<-1000
for (k in 1:n.discovery){
#zeroes trick for MLE~cond prob
tau[k]<- pow(SE[k], -2)
L[k]<- dnorm(MLE[k],mu.p53, tau[k])/(pnorm(-q*SE, mu.p53, tau[k]) +
1-pnorm(q*SE, mu.p53, tau[k]))
phi[k]<- -log(L[k])+C
zeroes[k]~dpois(phi[k])
}

mu.p53 ~ dnorm(0,.1)

}

cond.likelihood.re.model = function() {
for (j in 1:J) {
CaseCon[j] ~ dbern(theta[j])
logit(theta[j]) <-  beta.p53[site[j]]  }

C<-1000
for (k in 1:n.discovery){
#zeroes trick for MLE~cond prob
tau[k]<- pow(SE[k], -2)
L[k]<- dnorm(MLE[k],beta.p53[discovery.sites[k]], tau[k])/
(pnorm(-q*SE[k], beta.p53[discovery.sites[k]], tau[k]) +
1-pnorm(q*SE[k], beta.p53[discovery.sites[k]], tau[k]))
phi[k]<- -log(L[k])+C
zeroes[k]~dpois(phi[k])
}
for (i in 1:n.sites){ #total sites
beta.p53[i] ~dnorm(mu.p53,phi.p53)
}

mu.p53 ~ dnorm(0,.1)
# phi.p53 ~ dgamma(1, .05)
# sigma.p53 <- pow(phi.p53, -.5)
phi.p53 <- pow(sigma.p53, -2)
sigma.p53 ~ dt(0,1,1)%_%T(0,)
}

ones.cauchy.model= function() {
for (j in 1:J) {
CaseCon[j] ~ dbern(theta[j])
logit(theta[j]) <-  beta.p53[site[j]]  }

for (l in 1:n.sites) {
beta.p53.1[l] ~ dnorm(mu.p53, phi.p53)
beta.p53[l] <- beta.p53.1[l]*(mu.p53.notzero)
}
#ones trick
C<-1e6
epsilon<-0.01
tau<- pow(epsilon,-2)
mu.p53.notzero<- step(temp)
temp<-dt(mu.p53,0,1, 1)-dnorm(mu.p53,0,tau)
#cauchy is same as t with df=1
L<-(dnorm(mu.p53,0,tau)*(1-pind))+(dt(mu.p53,0,1,1)*pind)
p <- L/ C
one~ dbern(p)
mu.p53 ~ dunif(-10,10)
# phi.p53 ~ dgamma(1, .05)
# sigma.p53 <- pow(phi.p53, -.5)
phi.p53 <- pow(sigma.p53, -2)
sigma.p53 ~ dt(0,1,1)%_%T(0,)
pind ~ dbeta(.5,.5)
}

ones.normal.model= function() {
for (j in 1:J) {
CaseCon[j] ~ dbern(theta[j])
logit(theta[j]) <-  beta.p53[site[j]]  }

for (l in 1:n.sites) {
beta.p53.1[l] ~ dnorm(mu.p53, phi.p53)
beta.p53[l] <- beta.p53.1[l]*(mu.p53.notzero)
}
#ones trick
C<-1e6
epsilon<-0.01
tau<- pow(epsilon,-2)
mu.p53.notzero<- step(temp)
temp<-dnorm(mu.p53,0,1)-dnorm(mu.p53,0,tau)
L<-(dnorm(mu.p53,0,tau)*(1-pind))+(dnorm(mu.p53,0,1)*pind)
p <- L/ C
one ~ dbern(p)
mu.p53 ~ dunif(-10,10)
# phi.p53 ~ dgamma(1, .05)
# sigma.p53 <- pow(phi.p53, -.5)
phi.p53 <- pow(sigma.p53, -2)
sigma.p53 ~ dt(0,1,1)%_%T(0,)
pind ~ dbeta(.5,.5)
}

latent.normal.model= function() {
for (j in 1:J) {
CaseCon[j] ~ dbern(theta[j])
logit(theta[j]) <-  beta.p53[site[j]]  }

for (l in 1:n.sites) {
beta.p53.1[l] ~ dnorm(mu.p53, phi.p53)
beta.p53[l] <- beta.p53.1[l]*(mu.p53.notzero)
}

mu.p53<- mu1.p53*mu.p53.notzero
mu1.p53 ~ dnorm(0,1)
# phi.p53 ~ dgamma(1, .05)
# sigma.p53 <- pow(phi.p53, -.5)
phi.p53 <- pow(sigma.p53, -2)
sigma.p53 ~ dt(0,1,1)%_%T(0,)
mu.p53.notzero~dbern(pind)
pind ~ dbeta(.5,.5)
}

latent.cauchy.model= function() {
for (j in 1:J) {
CaseCon[j] ~ dbern(theta[j])
logit(theta[j]) <-  beta.p53[site[j]]  }

for (l in 1:n.sites) {
beta.p53.1[l] ~ dnorm(mu.p53, phi.p53)
beta.p53[l] <- beta.p53.1[l]*(mu.p53.notzero)
}

mu.p53<- mu1.p53*mu.p53.notzero
mu1.p53 ~ dt(0,1, 1)
# phi.p53 ~ dgamma(1, .05)
# sigma.p53 <- pow(phi.p53, -.5) #half cauchy, uniform to 1 or to 5, take a guess on sd
phi.p53 <- pow(sigma.p53, -2)
sigma.p53 ~ dt(0,1,1)%_%T(0,)
mu.p53.notzero~dbern(pind)
pind ~ dbeta(.5,.5)
}

fixed.model = function() {
for (j in 1:J) {
CaseCon[j] ~ dbern(theta[j])
logit(theta[j]) <-  mu.p53
}
mu.p53<- mu1.p53*mu.p53.notzero
mu1.p53 ~ dnorm(0,.1)
mu.p53.notzero~dbern(pind)
pind ~ dbeta(.5,.5)
}

original.model = function() {
for (j in 1:J) {
CaseCon[j] ~ dbern(theta[j])
logit(theta[j]) <-  beta.p53[site[j]]  }

for (l in 1:n.sites) {
beta.p53[l] ~ dnorm(mu.p53, phi.p53)
}

mu.p53 ~ dnorm(0,1)
phi.p53 ~ dgamma(1, .05)
sigma.p53 <- pow(phi.p53, -.5)
}

bf.model = function() {
for (j in 1:J) {
CaseCon[j] ~ dbern(theta[j])
logit(theta[j]) <-  beta.p53[site[j]] }

for (l in 1:n.sites) {
beta.p53.1[l] ~ dnorm(mu.p53, phi.p53)
beta.p53[l] <- beta.p53.1[l]*(mu.p53.notzero)
}

mu.p53<- mu1.p53*mu.p53.notzero
mu1.p53 ~ dt(0,1,1)
phi.p53 <- pow(sigma.p53, -2)
sigma.p53 ~ dt(0,1,1)%_%T(0,)

mu.p53.notzero~dbern(pind)
prior.odds<- (1-prior.pind)/prior.pind
BF<- 1/(-exp(1)*p*log(p)) #eplogp is BF h0/h1
pind<-prior.odds*BF/(1+prior.odds*BF)
prior.pind ~ dbeta(.9,.9)
}

## 7.2 Joint Analysis of TP53 SNPs

Table 7.1: Joint Analysis Estimates of $$\mu_{p53}$$
estimate 95% CI
rs9894946n 0 -0.018 - 0.172
rs12951053n 0 0 - 0 , 0.097 - 0.252
rs1625895n 0 0 - 0
rs1042522n 0 0 - 0
rs2909430n 0 -0.049 - 0.295
rs8079544n 0 0 - 0
rs2287499n 0 0 - 0 , 0.18 - 0.281
rs2287498n 0.27 0.083 - 0.445
rs8073498n 0 0 - 0
rs2078486n 0.372 -0.119 - -0.006 , 0 - 0.8