Bayesian HW9

1

fishdat <- read.csv("~/Desktop/fish.csv")

y <- as.factor(fishdat$survive)
y <- as.numeric(y)-1
x <- as.numeric(fishdat$duration)

x_sd <- sd(x)
x_mean <- mean(x)

x <- scale(x)
dat_scale <- as.data.frame(cbind(y,x))
glm(y~x, family = binomial(link="logit"))

## 
## Call:  glm(formula = y ~ x, family = binomial(link = "logit"))
## 
## Coefficients:
## (Intercept)            x  
##      0.7962       1.3212  
## 
## Degrees of Freedom: 159 Total (i.e. Null);  158 Residual
## Null Deviance:       209.5 
## Residual Deviance: 164.7     AIC: 168.7

set.seed(100)
tau2 <- 10
S <- 10000
propsd <- c(1,1)
x <- cbind(1, x)
p <- ncol(x)
beta <- rnorm(p)/10
Xbeta <- x%*% beta
beta_keep <- matrix(NA, S, p)


for(s in 1:S){
  for( j in 1:p){
    beta_prop <- rnorm(1)/propsd[j] + beta[j]
    Xbetaprop <- Xbeta + x[,j] * (beta_prop - beta[j])
    kprop <- sum(dbinom(y, 1, exp(Xbetaprop)/(1 + exp(Xbetaprop)), 
                        log = TRUE))
    kold <- sum(dbinom(y, 1, exp(Xbeta)/(1+exp(Xbeta)), 
                       log = TRUE))
    logR <- ((kprop - beta_prop^2/(2*tau2)) 
             - (kold - beta[j]^2/(2*tau2)))
    
    if(log(runif(1)) < logR){
      beta[j] <- beta_prop
      Xbeta <- Xbetaprop
    }
  }
  beta_keep [s,] <- beta
  
}

## trace Plots

matplot(beta_keep, type = 'l', las = 1)

#burnin

beta_keep <- beta_keep[-c(1:(S/2)),]

# acceptance ratio
colMeans(beta_keep[-1,] != beta_keep[-nrow(beta_keep),])

## [1] 0.2330466 0.2632527

## odds ratios

OR <- (beta_keep)

ORsummary <-
  data.frame(names = c("Beta0", "Beta1"),
             Mean = colMeans(OR),
             StdDev = apply(OR, 2, sd, 0.025),
             Lower = apply(OR, 2, quantile, 0.025),
             Upper = apply(OR, 2, quantile, 0.975))
ORsummary

##   names      Mean    StdDev     Lower    Upper
## 1 Beta0 0.8089137 0.2084994 0.4067013 1.233099
## 2 Beta1 1.3469732 0.2325737 0.9288680 1.834795

2

tumors <- read.csv("~/Desktop/tumors.csv")

x <- tumors$strain
x <- ifelse(x=="A", 1, 0)
X <- cbind(rep(1, length(x)), x)
y <- tumors$tumors

n <- length(y) 
p <- dim(X)[2]

beta_prior <- rep(0, p) 
beta_prop <- rep (10^2, p) 

library(MASS)

var_prop <- var(log(y+1/2))*solve(t(X)%*%X)
S<- 10000
beta<- rep (0,p)
acs<- 0
BETA<-matrix (0, nrow=S, ncol=p)

for(s in 1:S){
  beta.p <- t(mvrnorm(1 , beta , var_prop))
  beta.p <- as.numeric(beta.p)
  
  lhr <- sum(dpois(y,exp(X%*%beta.p), log=T)) -
    sum(dpois(y,exp(X%*%beta), log=T)) +
    sum(dnorm(beta.p,beta_prior,beta_prop, log=T)) -
    sum(dnorm(beta,beta_prior,beta_prop, log=T))
  
  if(log(runif(1)) < lhr ){ 
    beta<- beta.p
    acs<- acs+1
  }
  BETA[s,]<- beta
}

RR <- exp(BETA[,2])
results <- data.frame(
  mean = c(mean(BETA[,1]), mean(BETA[,2]), mean(RR), mean(RR>1)),
  sd = c(sd(BETA[,1]), sd(BETA[,2]), sd(RR), sd(RR>1)),
  quantile_lower = c(quantile(BETA[,1], 0.025), quantile(BETA[,2], 0.025),
                     quantile(RR, 0.025), quantile(RR>1, 0.025)),
  quantile_upper = c(quantile(BETA[,1], 0.975), quantile(BETA[,2], 0.975),
                     quantile(RR, 0.975), quantile(RR>1, 0.975)) 
)
rownames(results) <- c("beta0","beta1","RR", "RR>1")
results

##            mean        sd quantile_lower quantile_upper
## beta0 2.1459608 0.1798986     1.94174303      2.3451785
## beta1 0.2904366 0.1358652     0.02494026      0.5536439
## RR    1.3494216 0.1846399     1.02525387      1.7395803
## RR>1  0.9849000 0.1219569     1.00000000      1.0000000

## 2c ##
##  Number of Tumors for each Strain
strainA <- exp((BETA[,1]) + (BETA[,2]))
strainB <- exp((BETA[,1]))
mean(strainA)

## [1] 11.56427

mean(strainB)

## [1] 8.650882

### Prob strain A has more tumors than strain B
sum(strainA>strainB)/S

## [1] 0.9849

## this is the same as RR>1