• Bayesian Multilevel Logistic Regression
    • build the model and run the Gibbs sampler
  • Normal hierarchichal model
library(tidyr)
library("Przewodnik")
## Loading required package: PogromcyDanych
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: SmarterPoland
## Loading required package: ggplot2
## Loading required package: httr
## Loading required package: htmltools
## Loading required package: PBImisc
# library(dplyr)
library(ggplot2)
# library(tidyverse)
library(ggfortify)
library(e1071)
library(caTools)
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:httr':
## 
##     progress
library(brms)     # bayesian regression
## Loading required package: Rcpp
## Loading 'brms' package (version 2.17.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attaching package: 'brms'
## The following object is masked from 'package:e1071':
## 
##     rwiener
## The following object is masked from 'package:PBImisc':
## 
##     kidney
## The following object is masked from 'package:stats':
## 
##     ar
library(emmeans)    # used to get predicted means per condition
library(modelr)     # used to get predicted means per condition
## 
## Attaching package: 'modelr'
## The following object is masked from 'package:PBImisc':
## 
##     heights
library(tidybayes)  # for accessing model posteriors
## 
## Attaching package: 'tidybayes'
## The following objects are masked from 'package:brms':
## 
##     dstudent_t, pstudent_t, qstudent_t, rstudent_t
library(bayestestR)
## 
## Attaching package: 'bayestestR'
## The following object is masked from 'package:tidybayes':
## 
##     hdi

#Exploratory Dat Analysis

Loading the BRCA data

data(brca)

head(brca)
##   time death  MDM2    TP53 TP63      TP73 DNAJB1    DNAJB2 DNAJB4 DNAJB5 DNAJB6
## 1  259     0  0.22  1.4185 5.69  1.895648  -0.34 -0.389773   0.51  -0.10   0.16
## 2 1321     0  0.06 -0.2115 4.94  1.055648   0.25  0.340227  -0.21  -0.64   0.37
## 3 1463     0 -0.11  0.0685 4.23  0.645648   0.08  0.140227  -0.78  -0.71   0.10
## 4  434     0  1.29 -0.0315 1.43  1.465648   0.30  0.510227  -0.70   0.08   0.33
## 5 1437     0  0.14 -0.4315 3.81 -2.474352  -0.28  0.250227   0.04   0.30   0.24
## 6  635     0 -0.53  0.0485 1.42  1.575648  -0.11 -0.019773   0.35   1.11   0.74
##   DNAJB9   DNAJB11 DNAJB12     subtype p53mut
## 1  -0.11 -0.565677    0.07 BRCA.Normal      0
## 2  -0.18 -0.615677   -0.13   BRCA.LumA      0
## 3  -0.29 -0.505677    0.37   BRCA.LumA      0
## 4  -0.30 -0.795677    0.44   BRCA.LumA      0
## 5  -0.06 -0.005677    0.33   BRCA.LumA      0
## 6  -0.05  0.244323   -0.12   BRCA.LumB      1
summary(brca)
##       time            death             MDM2              TP53        
##  Min.   :   1.0   Min.   :0.0000   Min.   :-1.4700   Min.   :-3.0915  
##  1st Qu.: 477.5   1st Qu.:0.0000   1st Qu.:-0.3300   1st Qu.:-0.5415  
##  Median : 867.0   Median :0.0000   Median : 0.0100   Median :-0.0615  
##  Mean   :1224.3   Mean   :0.1813   Mean   : 0.1035   Mean   :-0.1662  
##  3rd Qu.:1627.0   3rd Qu.:0.0000   3rd Qu.: 0.4100   3rd Qu.: 0.2985  
##  Max.   :7126.0   Max.   :1.0000   Max.   : 4.6000   Max.   : 1.9185  
##                                    NA's   :50        NA's   :50       
##       TP63              TP73             DNAJB1             DNAJB2        
##  Min.   :-3.4300   Min.   :-4.0544   Min.   :-1.67000   Min.   :-2.50977  
##  1st Qu.:-0.1825   1st Qu.:-1.5643   1st Qu.:-0.33000   1st Qu.:-0.40977  
##  Median : 1.7350   Median :-0.5644   Median :-0.05000   Median :-0.05977  
##  Mean   : 1.4734   Mean   :-0.7029   Mean   :-0.01394   Mean   :-0.07627  
##  3rd Qu.: 3.3225   3rd Qu.: 0.2356   3rd Qu.: 0.29000   3rd Qu.: 0.29023  
##  Max.   : 6.6400   Max.   : 3.5957   Max.   : 2.28000   Max.   : 2.07023  
##  NA's   :67        NA's   :82        NA's   :50         NA's   :50        
##      DNAJB4             DNAJB5            DNAJB6            DNAJB9        
##  Min.   :-2.31000   Min.   :-2.5400   Min.   :-1.1700   Min.   :-2.01000  
##  1st Qu.:-0.51000   1st Qu.:-0.8200   1st Qu.: 0.0500   1st Qu.:-0.51000  
##  Median :-0.04000   Median :-0.3500   Median : 0.3400   Median :-0.13000  
##  Mean   :-0.04771   Mean   :-0.3739   Mean   : 0.3587   Mean   :-0.07747  
##  3rd Qu.: 0.37000   3rd Qu.: 0.1200   3rd Qu.: 0.6500   3rd Qu.: 0.28000  
##  Max.   : 1.84000   Max.   : 2.1600   Max.   : 2.1800   Max.   : 2.54000  
##  NA's   :50         NA's   :50        NA's   :50        NA's   :50        
##     DNAJB11           DNAJB12                subtype        p53mut      
##  Min.   :-1.3957   Min.   :-1.90000   BRCA.Basal :109   Min.   :0.0000  
##  1st Qu.:-0.6257   1st Qu.:-0.17000   BRCA.Her2  : 43   1st Qu.:0.0000  
##  Median :-0.2557   Median : 0.08000   BRCA.LumA  :365   Median :0.0000  
##  Mean   :-0.1785   Mean   : 0.05896   BRCA.LumB  :156   Mean   :0.3304  
##  3rd Qu.: 0.1643   3rd Qu.: 0.31000   BRCA.Normal: 22   3rd Qu.:1.0000  
##  Max.   : 2.0443   Max.   : 1.99000                     Max.   :1.0000  
##  NA's   :50        NA's   :50                           NA's   :126
library(dplyr)

brca.genes = select(brca, -c(subtype, time, death, p53mut))

plot means

boxplot(brca.genes, las = 2, col = terrain.colors(12),
        main = "12 gene expression profiles from BRCA dataset",
        xlab = "", 
        ylab = "Gene expression")

mtext("Genes", side=1, line = -1)

# dev.off()
brc.genes_plot = ggplot(brca.genes, aes(x =as.factor(brca.genes), y = expression)) +
  geom_boxplot()
unique(brca$subtype)
## [1] BRCA.Normal BRCA.LumA   BRCA.LumB   BRCA.Basal  BRCA.Her2  
## Levels: BRCA.Basal BRCA.Her2 BRCA.LumA BRCA.LumB BRCA.Normal

number of normal

n.norm = filter(brca, subtype == "BRCA.Normal")
nrow(n.norm)
## [1] 22

number of lumunal A

n.lumA = filter(brca, subtype == "BRCA.LumA")
nrow(n.lumA)
## [1] 365

number of luminal B

n.lumB = filter(brca, subtype == "BRCA.LumB")
nrow(n.lumB)
## [1] 156

number of her2

n.her2 = filter(brca, subtype == "BRCA.Her2")
nrow(n.her2)
## [1] 43

number of basal

n.basal = filter(brca, subtype == "BRCA.Basal")
nrow(n.basal)
## [1] 109
colnames(brca)
##  [1] "time"    "death"   "MDM2"    "TP53"    "TP63"    "TP73"    "DNAJB1" 
##  [8] "DNAJB2"  "DNAJB4"  "DNAJB5"  "DNAJB6"  "DNAJB9"  "DNAJB11" "DNAJB12"
## [15] "subtype" "p53mut"
pca2 = prcomp(drop_na(brca[, c(3:14)]))

summary(pca2)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.2335 1.3571 0.82621 0.76011 0.72231 0.64642 0.55270
## Proportion of Variance 0.4807 0.1775 0.06578 0.05567 0.05027 0.04026 0.02944
## Cumulative Proportion  0.4807 0.6582 0.72393 0.77961 0.82988 0.87014 0.89958
##                            PC8     PC9    PC10   PC11    PC12
## Standard deviation     0.53215 0.47058 0.45514 0.4225 0.38968
## Proportion of Variance 0.02729 0.02134 0.01996 0.0172 0.01463
## Cumulative Proportion  0.92687 0.94820 0.96816 0.9854 1.00000
str(pca2)
## List of 5
##  $ sdev    : num [1:12] 2.234 1.357 0.826 0.76 0.722 ...
##  $ rotation: num [1:12, 1:12] -0.0122 -0.0081 -0.9905 -0.0661 0.0119 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:12] "MDM2" "TP53" "TP63" "TP73" ...
##   .. ..$ : chr [1:12] "PC1" "PC2" "PC3" "PC4" ...
##  $ center  : Named num [1:12] 0.1322 -0.1357 1.5898 -0.6946 -0.0288 ...
##   ..- attr(*, "names")= chr [1:12] "MDM2" "TP53" "TP63" "TP73" ...
##  $ scale   : logi FALSE
##  $ x       : num [1:598, 1:12] -4.3373 -3.4081 -2.6446 -0.0103 -2.1388 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:12] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "prcomp"

PCA plot of all subtypes

pca2.plot = autoplot(pca2, data = drop_na(brca[, c(3:15)]), colour = 'subtype')

pca2.plot

colnames(drop_na(brca[, c(3:15)]))
##  [1] "MDM2"    "TP53"    "TP63"    "TP73"    "DNAJB1"  "DNAJB2"  "DNAJB4" 
##  [8] "DNAJB5"  "DNAJB6"  "DNAJB9"  "DNAJB11" "DNAJB12" "subtype"

PCA subset of two groups: Luminal A vs Luminal B

brca.lumAB = filter(brca, subtype == "BRCA.LumA" | subtype == "BRCA.LumB")

nrow(brca.lumAB)
## [1] 521
head(brca.lumAB)
##   time death  MDM2    TP53 TP63      TP73 DNAJB1    DNAJB2 DNAJB4 DNAJB5 DNAJB6
## 1 1321     0  0.06 -0.2115 4.94  1.055648   0.25  0.340227  -0.21  -0.64   0.37
## 2 1463     0 -0.11  0.0685 4.23  0.645648   0.08  0.140227  -0.78  -0.71   0.10
## 3  434     0  1.29 -0.0315 1.43  1.465648   0.30  0.510227  -0.70   0.08   0.33
## 4 1437     0  0.14 -0.4315 3.81 -2.474352  -0.28  0.250227   0.04   0.30   0.24
## 5  635     0 -0.53  0.0485 1.42  1.575648  -0.11 -0.019773   0.35   1.11   0.74
## 6  416     0  3.17 -0.6015 3.54 -0.504352  -0.20 -0.209773   0.08  -0.59   0.53
##   DNAJB9   DNAJB11 DNAJB12   subtype p53mut
## 1  -0.18 -0.615677   -0.13 BRCA.LumA      0
## 2  -0.29 -0.505677    0.37 BRCA.LumA      0
## 3  -0.30 -0.795677    0.44 BRCA.LumA      0
## 4  -0.06 -0.005677    0.33 BRCA.LumA      0
## 5  -0.05  0.244323   -0.12 BRCA.LumB      1
## 6   0.53  0.064323   -0.39 BRCA.LumA      0
pca3 = prcomp(drop_na(brca.lumAB[, c(3:14)]))

summary(pca3)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.1920 1.2318 0.79217 0.68256 0.63148 0.60793 0.53417
## Proportion of Variance 0.5159 0.1629 0.06738 0.05002 0.04282 0.03968 0.03064
## Cumulative Proportion  0.5159 0.6789 0.74625 0.79627 0.83909 0.87877 0.90941
##                            PC8    PC9    PC10    PC11    PC12
## Standard deviation     0.46753 0.4454 0.41183 0.37592 0.34026
## Proportion of Variance 0.02347 0.0213 0.01821 0.01517 0.01243
## Cumulative Proportion  0.93288 0.9542 0.97240 0.98757 1.00000
str(pca3)
## List of 5
##  $ sdev    : num [1:12] 2.192 1.232 0.792 0.683 0.631 ...
##  $ rotation: num [1:12, 1:12] 0.00951 0.00996 -0.99111 -0.00487 0.01577 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:12] "MDM2" "TP53" "TP63" "TP73" ...
##   .. ..$ : chr [1:12] "PC1" "PC2" "PC3" "PC4" ...
##  $ center  : Named num [1:12] 0.2526 -0.1415 1.6752 -0.5133 -0.0431 ...
##   ..- attr(*, "names")= chr [1:12] "MDM2" "TP53" "TP63" "TP73" ...
##  $ scale   : logi FALSE
##  $ x       : num [1:471, 1:12] -3.2127 -2.457 0.2592 -2.1713 0.0867 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:12] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "prcomp"
pca3.plot = autoplot(pca3, data = drop_na(brca.lumAB[, c(3:15)]), colour = 'subtype')

pca3.plot

Bayesian Multilevel Logistic Regression

build the model and run the Gibbs sampler

logit_model = brm(subtype ~., data = brca.lumAB, family = bernoulli, prior=set_prior("normal(0,2)", class = "b"))
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include   -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## 
## SAMPLING FOR MODEL 'bde997d76b8897c7d3055547ebdfb0af' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 2.12206 seconds (Warm-up)
## Chain 1:                0.873771 seconds (Sampling)
## Chain 1:                2.99583 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'bde997d76b8897c7d3055547ebdfb0af' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 2.16999 seconds (Warm-up)
## Chain 2:                0.705839 seconds (Sampling)
## Chain 2:                2.87583 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'bde997d76b8897c7d3055547ebdfb0af' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.6e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 2.18978 seconds (Warm-up)
## Chain 3:                0.862753 seconds (Sampling)
## Chain 3:                3.05253 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'bde997d76b8897c7d3055547ebdfb0af' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.5e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2.23668 seconds (Warm-up)
## Chain 4:                0.76633 seconds (Sampling)
## Chain 4:                3.00301 seconds (Total)
## Chain 4:

Summarize the model

summary(logit_model)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: subtype ~ time + death + MDM2 + TP53 + TP63 + TP73 + DNAJB1 + DNAJB2 + DNAJB4 + DNAJB5 + DNAJB6 + DNAJB9 + DNAJB11 + DNAJB12 + p53mut 
##    Data: brca.lumAB (Number of observations: 409) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.20      0.34    -0.88     0.46 1.00     3459     2934
## time         -0.00      0.00    -0.00     0.00 1.00     3863     3087
## death         0.03      0.44    -0.83     0.90 1.00     3327     3065
## MDM2          0.18      0.26    -0.31     0.68 1.00     3148     3201
## TP53          0.20      0.23    -0.27     0.64 1.00     3080     2653
## TP63         -0.65      0.08    -0.82    -0.48 1.00     2976     3065
## TP73          0.35      0.13     0.09     0.62 1.00     3140     3036
## DNAJB1        0.40      0.35    -0.30     1.10 1.00     3505     2593
## DNAJB2       -1.32      0.29    -1.89    -0.76 1.00     3352     2593
## DNAJB4       -0.22      0.26    -0.74     0.28 1.00     3288     2799
## DNAJB5       -0.87      0.24    -1.36    -0.41 1.00     2923     2833
## DNAJB6        0.32      0.33    -0.32     0.96 1.00     3459     3011
## DNAJB9       -0.54      0.28    -1.08    -0.00 1.00     3136     3088
## DNAJB11       1.54      0.35     0.89     2.25 1.00     2640     2868
## DNAJB12       0.00      0.39    -0.76     0.74 1.00     3487     3028
## p53mut        1.24      0.40     0.47     2.02 1.00     2989     2863
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Describe the posterior

describe_posterior(logit_model, ci = .95, rope_ci = 0.95, test = c('pd', 'p_map', 'rope', 'bf'))
## Sampling priors, please wait...
## Summary of Posterior Distribution
## 
## Parameter   |    Median |         95% CI | p (MAP) |     pd |          ROPE | % in ROPE |  Rhat |     ESS |      BF
## -------------------------------------------------------------------------------------------------------------------
## (Intercept) |     -0.20 | [-0.89,  0.44] |  0.866  | 71.40% | [-0.18, 0.18] |    36.81% | 1.000 | 3431.00 | < 0.001
## time        | -1.89e-04 | [ 0.00,  0.00] |  0.463  | 89.92% | [-0.18, 0.18] |      100% | 1.001 | 3826.00 | < 0.001
## death       |      0.03 | [-0.84,  0.89] |  0.994  | 52.80% | [-0.18, 0.18] |    32.54% | 0.999 | 3309.00 |   0.221
## MDM2        |      0.18 | [-0.33,  0.66] |  0.756  | 75.83% | [-0.18, 0.18] |    44.09% | 1.000 | 3146.00 |   0.161
## TP53        |      0.21 | [-0.27,  0.63] |  0.642  | 81.50% | [-0.18, 0.18] |    42.07% | 1.000 | 3004.00 |   0.174
## TP63        |     -0.65 | [-0.82, -0.48] |  < .001 |   100% | [-0.18, 0.18] |        0% | 0.999 | 2943.00 |  > 1000
## TP73        |      0.35 | [ 0.08,  0.61] |  0.039  | 99.60% | [-0.18, 0.18] |     8.84% | 1.000 | 3146.00 |    1.75
## DNAJB1      |      0.40 | [-0.30,  1.09] |  0.597  | 86.55% | [-0.18, 0.18] |    23.86% | 1.000 | 3482.00 |   0.314
## DNAJB2      |     -1.32 | [-1.89, -0.76] |  < .001 |   100% | [-0.18, 0.18] |        0% | 1.000 | 3322.00 |  513.14
## DNAJB4      |     -0.22 | [-0.71,  0.29] |  0.686  | 80.65% | [-0.18, 0.18] |    40.33% | 1.000 | 3261.00 |   0.180
## DNAJB5      |     -0.86 | [-1.32, -0.39] |  < .001 |   100% | [-0.18, 0.18] |        0% | 1.000 | 2905.00 |   99.10
## DNAJB6      |      0.31 | [-0.35,  0.92] |  0.680  | 83.50% | [-0.18, 0.18] |    29.44% | 1.000 | 3443.00 |   0.248
## DNAJB9      |     -0.54 | [-1.10, -0.04] |  0.128  | 97.50% | [-0.18, 0.18] |     7.00% | 1.000 | 3119.00 |    1.05
## DNAJB11     |      1.53 | [ 0.88,  2.24] |  < .001 |   100% | [-0.18, 0.18] |        0% | 1.000 | 2614.00 |  > 1000
## DNAJB12     |  6.78e-03 | [-0.77,  0.74] |  0.997  | 50.65% | [-0.18, 0.18] |    38.28% | 1.000 | 3452.00 |   0.199
## p53mut      |      1.24 | [ 0.47,  2.02] |  < .001 |   100% | [-0.18, 0.18] |        0% | 1.000 | 2961.00 |   67.45

PLot the posterior estimtes of the model parameters

plot(logit_model)

summary(logit_model)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: subtype ~ time + death + MDM2 + TP53 + TP63 + TP73 + DNAJB1 + DNAJB2 + DNAJB4 + DNAJB5 + DNAJB6 + DNAJB9 + DNAJB11 + DNAJB12 + p53mut 
##    Data: brca.lumAB (Number of observations: 409) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.20      0.34    -0.88     0.46 1.00     3459     2934
## time         -0.00      0.00    -0.00     0.00 1.00     3863     3087
## death         0.03      0.44    -0.83     0.90 1.00     3327     3065
## MDM2          0.18      0.26    -0.31     0.68 1.00     3148     3201
## TP53          0.20      0.23    -0.27     0.64 1.00     3080     2653
## TP63         -0.65      0.08    -0.82    -0.48 1.00     2976     3065
## TP73          0.35      0.13     0.09     0.62 1.00     3140     3036
## DNAJB1        0.40      0.35    -0.30     1.10 1.00     3505     2593
## DNAJB2       -1.32      0.29    -1.89    -0.76 1.00     3352     2593
## DNAJB4       -0.22      0.26    -0.74     0.28 1.00     3288     2799
## DNAJB5       -0.87      0.24    -1.36    -0.41 1.00     2923     2833
## DNAJB6        0.32      0.33    -0.32     0.96 1.00     3459     3011
## DNAJB9       -0.54      0.28    -1.08    -0.00 1.00     3136     3088
## DNAJB11       1.54      0.35     0.89     2.25 1.00     2640     2868
## DNAJB12       0.00      0.39    -0.76     0.74 1.00     3487     3028
## p53mut        1.24      0.40     0.47     2.02 1.00     2989     2863
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Store the posterior samples

post = posterior_samples(logit_model)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.

Calculate point estimates of the posterior samples

mean(post$b_TP53) ## high in luminal A
## [1] 0.2026549
mean(post$b_MDM2) #general luminal specific marker
## [1] 0.1839254
print(c(mean(post$b_TP53), mean(post$b_MDM2)))
## [1] 0.2026549 0.1839254

(Optional) Perform an evaluation of the group means

# install.packages("magrittr")
# install.packages("dplyr", dependencies = T) 
# library(magrittr)
# library(tidyverse)
# library(dplyr)

# invlogit = function(x){
#   exp(x) / (1 + exp(x))
# }

# library(tidyverse)
# library(dplyr)

# post %>% mutate(p_0.3_0.2 = invlogit(b_Intercept + b_TP53*0.3 + b_MDM2*0.2), p_0.05_0.08 = invlogit(b_Intercept + b_TP53*0.05 + b_MDM2*0.08)) = post


# post$p_0.3_0.2 = invlogit(b_Intercept + b_TP53*0.3 + b_MDM2*0.2)
# post$p_0.05_0.08 = invlogit(b_Intercept + b_TP53*0.05 + b_MDM2*0.08)
pp_check(logit_model, ndraws = 1000,
         resp = subtype)

Normal hierarchichal model

brca.lumAB.genes  = filter(brca.lumAB[3:14])

par(mfrow=c(1,1))
boxplot(select(brca, -c(subtype, time, death, p53mut)))

Calculating necessary statistics:

##Calculating necessary statistics:
y = drop_na(brca.lumAB.genes)

m = length(y[1,])
n = rep(NA,m)
means = rep(NA,m)
for (i in 1:m){
  n[i] = length(y[,i])
  means[i] = mean(y[,i])
}
ntot = sum(n)


## true sigmasq
truesigsq = 8 

Sampling Parameters for Normal Hierarchical Model

#####################################################
# Sampling Parameters for Normal Hierarchical Model #
#####################################################

## finding right grid for tausq
tausq.grid = ppoints(1000)*20

tausq.logpostfunc = function(tausq){
    Vmu0 = 1/sum(1/(tausq + truesigsq/n))
    mu0hat = sum(means/(tausq + truesigsq/n))*Vmu0
    out = -0.5*log(tausq)+0.5*log(Vmu0)
    for (group in 1:m){
        out = out - 0.5*log(tausq + truesigsq/n[group])
    }
    for (group in 1:m){
        out = out - 0.5*((means[group]-mu0hat)^2)/(tausq + truesigsq/n[group])
    }
    out
}
tausq.logpost = rep(NA,1000)
for (i in 1:1000){
    tausq.logpost[i] = tausq.logpostfunc(tausq.grid[i])
}
tausq.post = exp(tausq.logpost-max(tausq.logpost))
tausq.post = tausq.post/sum(tausq.post)

par(mfrow=c(1,1))
plot(tausq.grid,tausq.post,type="l")

numsamp = 1000
tausq.samp = rep(NA,numsamp)
mu0.samp = rep(NA,numsamp)
mu.samp = matrix(NA,nrow=numsamp,ncol=m)
for (i in 1:numsamp){
    # sampling tausq from grid of values
    curtausq = sample(tausq.grid,size=1,prob=tausq.post)
    # sampling mu0 given curtausq   
    Vmu0 = 1/sum(1/(curtausq + truesigsq/n))
    mu0hat = sum(means/(curtausq + truesigsq/n))*Vmu0
    curmu0 = rnorm(1,mean=mu0hat,sd=sqrt(Vmu0))
    # sampling group means given curtausq and curmu0
    curmu = rep(NA,m)
    for (j in 1:m){
        curvar = 1/(n[j]/truesigsq + 1/curtausq)
        curmean = (means[j]*n[j]/truesigsq + curmu0/curtausq)*curvar
        curmu[j] = rnorm(1,mean=curmean,sd=sqrt(curvar))
    }
    tausq.samp[i] = curtausq
    mu0.samp[i] = curmu0
    mu.samp[i,] = curmu
    print (i)
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32
## [1] 33
## [1] 34
## [1] 35
## [1] 36
## [1] 37
## [1] 38
## [1] 39
## [1] 40
## [1] 41
## [1] 42
## [1] 43
## [1] 44
## [1] 45
## [1] 46
## [1] 47
## [1] 48
## [1] 49
## [1] 50
## [1] 51
## [1] 52
## [1] 53
## [1] 54
## [1] 55
## [1] 56
## [1] 57
## [1] 58
## [1] 59
## [1] 60
## [1] 61
## [1] 62
## [1] 63
## [1] 64
## [1] 65
## [1] 66
## [1] 67
## [1] 68
## [1] 69
## [1] 70
## [1] 71
## [1] 72
## [1] 73
## [1] 74
## [1] 75
## [1] 76
## [1] 77
## [1] 78
## [1] 79
## [1] 80
## [1] 81
## [1] 82
## [1] 83
## [1] 84
## [1] 85
## [1] 86
## [1] 87
## [1] 88
## [1] 89
## [1] 90
## [1] 91
## [1] 92
## [1] 93
## [1] 94
## [1] 95
## [1] 96
## [1] 97
## [1] 98
## [1] 99
## [1] 100
## [1] 101
## [1] 102
## [1] 103
## [1] 104
## [1] 105
## [1] 106
## [1] 107
## [1] 108
## [1] 109
## [1] 110
## [1] 111
## [1] 112
## [1] 113
## [1] 114
## [1] 115
## [1] 116
## [1] 117
## [1] 118
## [1] 119
## [1] 120
## [1] 121
## [1] 122
## [1] 123
## [1] 124
## [1] 125
## [1] 126
## [1] 127
## [1] 128
## [1] 129
## [1] 130
## [1] 131
## [1] 132
## [1] 133
## [1] 134
## [1] 135
## [1] 136
## [1] 137
## [1] 138
## [1] 139
## [1] 140
## [1] 141
## [1] 142
## [1] 143
## [1] 144
## [1] 145
## [1] 146
## [1] 147
## [1] 148
## [1] 149
## [1] 150
## [1] 151
## [1] 152
## [1] 153
## [1] 154
## [1] 155
## [1] 156
## [1] 157
## [1] 158
## [1] 159
## [1] 160
## [1] 161
## [1] 162
## [1] 163
## [1] 164
## [1] 165
## [1] 166
## [1] 167
## [1] 168
## [1] 169
## [1] 170
## [1] 171
## [1] 172
## [1] 173
## [1] 174
## [1] 175
## [1] 176
## [1] 177
## [1] 178
## [1] 179
## [1] 180
## [1] 181
## [1] 182
## [1] 183
## [1] 184
## [1] 185
## [1] 186
## [1] 187
## [1] 188
## [1] 189
## [1] 190
## [1] 191
## [1] 192
## [1] 193
## [1] 194
## [1] 195
## [1] 196
## [1] 197
## [1] 198
## [1] 199
## [1] 200
## [1] 201
## [1] 202
## [1] 203
## [1] 204
## [1] 205
## [1] 206
## [1] 207
## [1] 208
## [1] 209
## [1] 210
## [1] 211
## [1] 212
## [1] 213
## [1] 214
## [1] 215
## [1] 216
## [1] 217
## [1] 218
## [1] 219
## [1] 220
## [1] 221
## [1] 222
## [1] 223
## [1] 224
## [1] 225
## [1] 226
## [1] 227
## [1] 228
## [1] 229
## [1] 230
## [1] 231
## [1] 232
## [1] 233
## [1] 234
## [1] 235
## [1] 236
## [1] 237
## [1] 238
## [1] 239
## [1] 240
## [1] 241
## [1] 242
## [1] 243
## [1] 244
## [1] 245
## [1] 246
## [1] 247
## [1] 248
## [1] 249
## [1] 250
## [1] 251
## [1] 252
## [1] 253
## [1] 254
## [1] 255
## [1] 256
## [1] 257
## [1] 258
## [1] 259
## [1] 260
## [1] 261
## [1] 262
## [1] 263
## [1] 264
## [1] 265
## [1] 266
## [1] 267
## [1] 268
## [1] 269
## [1] 270
## [1] 271
## [1] 272
## [1] 273
## [1] 274
## [1] 275
## [1] 276
## [1] 277
## [1] 278
## [1] 279
## [1] 280
## [1] 281
## [1] 282
## [1] 283
## [1] 284
## [1] 285
## [1] 286
## [1] 287
## [1] 288
## [1] 289
## [1] 290
## [1] 291
## [1] 292
## [1] 293
## [1] 294
## [1] 295
## [1] 296
## [1] 297
## [1] 298
## [1] 299
## [1] 300
## [1] 301
## [1] 302
## [1] 303
## [1] 304
## [1] 305
## [1] 306
## [1] 307
## [1] 308
## [1] 309
## [1] 310
## [1] 311
## [1] 312
## [1] 313
## [1] 314
## [1] 315
## [1] 316
## [1] 317
## [1] 318
## [1] 319
## [1] 320
## [1] 321
## [1] 322
## [1] 323
## [1] 324
## [1] 325
## [1] 326
## [1] 327
## [1] 328
## [1] 329
## [1] 330
## [1] 331
## [1] 332
## [1] 333
## [1] 334
## [1] 335
## [1] 336
## [1] 337
## [1] 338
## [1] 339
## [1] 340
## [1] 341
## [1] 342
## [1] 343
## [1] 344
## [1] 345
## [1] 346
## [1] 347
## [1] 348
## [1] 349
## [1] 350
## [1] 351
## [1] 352
## [1] 353
## [1] 354
## [1] 355
## [1] 356
## [1] 357
## [1] 358
## [1] 359
## [1] 360
## [1] 361
## [1] 362
## [1] 363
## [1] 364
## [1] 365
## [1] 366
## [1] 367
## [1] 368
## [1] 369
## [1] 370
## [1] 371
## [1] 372
## [1] 373
## [1] 374
## [1] 375
## [1] 376
## [1] 377
## [1] 378
## [1] 379
## [1] 380
## [1] 381
## [1] 382
## [1] 383
## [1] 384
## [1] 385
## [1] 386
## [1] 387
## [1] 388
## [1] 389
## [1] 390
## [1] 391
## [1] 392
## [1] 393
## [1] 394
## [1] 395
## [1] 396
## [1] 397
## [1] 398
## [1] 399
## [1] 400
## [1] 401
## [1] 402
## [1] 403
## [1] 404
## [1] 405
## [1] 406
## [1] 407
## [1] 408
## [1] 409
## [1] 410
## [1] 411
## [1] 412
## [1] 413
## [1] 414
## [1] 415
## [1] 416
## [1] 417
## [1] 418
## [1] 419
## [1] 420
## [1] 421
## [1] 422
## [1] 423
## [1] 424
## [1] 425
## [1] 426
## [1] 427
## [1] 428
## [1] 429
## [1] 430
## [1] 431
## [1] 432
## [1] 433
## [1] 434
## [1] 435
## [1] 436
## [1] 437
## [1] 438
## [1] 439
## [1] 440
## [1] 441
## [1] 442
## [1] 443
## [1] 444
## [1] 445
## [1] 446
## [1] 447
## [1] 448
## [1] 449
## [1] 450
## [1] 451
## [1] 452
## [1] 453
## [1] 454
## [1] 455
## [1] 456
## [1] 457
## [1] 458
## [1] 459
## [1] 460
## [1] 461
## [1] 462
## [1] 463
## [1] 464
## [1] 465
## [1] 466
## [1] 467
## [1] 468
## [1] 469
## [1] 470
## [1] 471
## [1] 472
## [1] 473
## [1] 474
## [1] 475
## [1] 476
## [1] 477
## [1] 478
## [1] 479
## [1] 480
## [1] 481
## [1] 482
## [1] 483
## [1] 484
## [1] 485
## [1] 486
## [1] 487
## [1] 488
## [1] 489
## [1] 490
## [1] 491
## [1] 492
## [1] 493
## [1] 494
## [1] 495
## [1] 496
## [1] 497
## [1] 498
## [1] 499
## [1] 500
## [1] 501
## [1] 502
## [1] 503
## [1] 504
## [1] 505
## [1] 506
## [1] 507
## [1] 508
## [1] 509
## [1] 510
## [1] 511
## [1] 512
## [1] 513
## [1] 514
## [1] 515
## [1] 516
## [1] 517
## [1] 518
## [1] 519
## [1] 520
## [1] 521
## [1] 522
## [1] 523
## [1] 524
## [1] 525
## [1] 526
## [1] 527
## [1] 528
## [1] 529
## [1] 530
## [1] 531
## [1] 532
## [1] 533
## [1] 534
## [1] 535
## [1] 536
## [1] 537
## [1] 538
## [1] 539
## [1] 540
## [1] 541
## [1] 542
## [1] 543
## [1] 544
## [1] 545
## [1] 546
## [1] 547
## [1] 548
## [1] 549
## [1] 550
## [1] 551
## [1] 552
## [1] 553
## [1] 554
## [1] 555
## [1] 556
## [1] 557
## [1] 558
## [1] 559
## [1] 560
## [1] 561
## [1] 562
## [1] 563
## [1] 564
## [1] 565
## [1] 566
## [1] 567
## [1] 568
## [1] 569
## [1] 570
## [1] 571
## [1] 572
## [1] 573
## [1] 574
## [1] 575
## [1] 576
## [1] 577
## [1] 578
## [1] 579
## [1] 580
## [1] 581
## [1] 582
## [1] 583
## [1] 584
## [1] 585
## [1] 586
## [1] 587
## [1] 588
## [1] 589
## [1] 590
## [1] 591
## [1] 592
## [1] 593
## [1] 594
## [1] 595
## [1] 596
## [1] 597
## [1] 598
## [1] 599
## [1] 600
## [1] 601
## [1] 602
## [1] 603
## [1] 604
## [1] 605
## [1] 606
## [1] 607
## [1] 608
## [1] 609
## [1] 610
## [1] 611
## [1] 612
## [1] 613
## [1] 614
## [1] 615
## [1] 616
## [1] 617
## [1] 618
## [1] 619
## [1] 620
## [1] 621
## [1] 622
## [1] 623
## [1] 624
## [1] 625
## [1] 626
## [1] 627
## [1] 628
## [1] 629
## [1] 630
## [1] 631
## [1] 632
## [1] 633
## [1] 634
## [1] 635
## [1] 636
## [1] 637
## [1] 638
## [1] 639
## [1] 640
## [1] 641
## [1] 642
## [1] 643
## [1] 644
## [1] 645
## [1] 646
## [1] 647
## [1] 648
## [1] 649
## [1] 650
## [1] 651
## [1] 652
## [1] 653
## [1] 654
## [1] 655
## [1] 656
## [1] 657
## [1] 658
## [1] 659
## [1] 660
## [1] 661
## [1] 662
## [1] 663
## [1] 664
## [1] 665
## [1] 666
## [1] 667
## [1] 668
## [1] 669
## [1] 670
## [1] 671
## [1] 672
## [1] 673
## [1] 674
## [1] 675
## [1] 676
## [1] 677
## [1] 678
## [1] 679
## [1] 680
## [1] 681
## [1] 682
## [1] 683
## [1] 684
## [1] 685
## [1] 686
## [1] 687
## [1] 688
## [1] 689
## [1] 690
## [1] 691
## [1] 692
## [1] 693
## [1] 694
## [1] 695
## [1] 696
## [1] 697
## [1] 698
## [1] 699
## [1] 700
## [1] 701
## [1] 702
## [1] 703
## [1] 704
## [1] 705
## [1] 706
## [1] 707
## [1] 708
## [1] 709
## [1] 710
## [1] 711
## [1] 712
## [1] 713
## [1] 714
## [1] 715
## [1] 716
## [1] 717
## [1] 718
## [1] 719
## [1] 720
## [1] 721
## [1] 722
## [1] 723
## [1] 724
## [1] 725
## [1] 726
## [1] 727
## [1] 728
## [1] 729
## [1] 730
## [1] 731
## [1] 732
## [1] 733
## [1] 734
## [1] 735
## [1] 736
## [1] 737
## [1] 738
## [1] 739
## [1] 740
## [1] 741
## [1] 742
## [1] 743
## [1] 744
## [1] 745
## [1] 746
## [1] 747
## [1] 748
## [1] 749
## [1] 750
## [1] 751
## [1] 752
## [1] 753
## [1] 754
## [1] 755
## [1] 756
## [1] 757
## [1] 758
## [1] 759
## [1] 760
## [1] 761
## [1] 762
## [1] 763
## [1] 764
## [1] 765
## [1] 766
## [1] 767
## [1] 768
## [1] 769
## [1] 770
## [1] 771
## [1] 772
## [1] 773
## [1] 774
## [1] 775
## [1] 776
## [1] 777
## [1] 778
## [1] 779
## [1] 780
## [1] 781
## [1] 782
## [1] 783
## [1] 784
## [1] 785
## [1] 786
## [1] 787
## [1] 788
## [1] 789
## [1] 790
## [1] 791
## [1] 792
## [1] 793
## [1] 794
## [1] 795
## [1] 796
## [1] 797
## [1] 798
## [1] 799
## [1] 800
## [1] 801
## [1] 802
## [1] 803
## [1] 804
## [1] 805
## [1] 806
## [1] 807
## [1] 808
## [1] 809
## [1] 810
## [1] 811
## [1] 812
## [1] 813
## [1] 814
## [1] 815
## [1] 816
## [1] 817
## [1] 818
## [1] 819
## [1] 820
## [1] 821
## [1] 822
## [1] 823
## [1] 824
## [1] 825
## [1] 826
## [1] 827
## [1] 828
## [1] 829
## [1] 830
## [1] 831
## [1] 832
## [1] 833
## [1] 834
## [1] 835
## [1] 836
## [1] 837
## [1] 838
## [1] 839
## [1] 840
## [1] 841
## [1] 842
## [1] 843
## [1] 844
## [1] 845
## [1] 846
## [1] 847
## [1] 848
## [1] 849
## [1] 850
## [1] 851
## [1] 852
## [1] 853
## [1] 854
## [1] 855
## [1] 856
## [1] 857
## [1] 858
## [1] 859
## [1] 860
## [1] 861
## [1] 862
## [1] 863
## [1] 864
## [1] 865
## [1] 866
## [1] 867
## [1] 868
## [1] 869
## [1] 870
## [1] 871
## [1] 872
## [1] 873
## [1] 874
## [1] 875
## [1] 876
## [1] 877
## [1] 878
## [1] 879
## [1] 880
## [1] 881
## [1] 882
## [1] 883
## [1] 884
## [1] 885
## [1] 886
## [1] 887
## [1] 888
## [1] 889
## [1] 890
## [1] 891
## [1] 892
## [1] 893
## [1] 894
## [1] 895
## [1] 896
## [1] 897
## [1] 898
## [1] 899
## [1] 900
## [1] 901
## [1] 902
## [1] 903
## [1] 904
## [1] 905
## [1] 906
## [1] 907
## [1] 908
## [1] 909
## [1] 910
## [1] 911
## [1] 912
## [1] 913
## [1] 914
## [1] 915
## [1] 916
## [1] 917
## [1] 918
## [1] 919
## [1] 920
## [1] 921
## [1] 922
## [1] 923
## [1] 924
## [1] 925
## [1] 926
## [1] 927
## [1] 928
## [1] 929
## [1] 930
## [1] 931
## [1] 932
## [1] 933
## [1] 934
## [1] 935
## [1] 936
## [1] 937
## [1] 938
## [1] 939
## [1] 940
## [1] 941
## [1] 942
## [1] 943
## [1] 944
## [1] 945
## [1] 946
## [1] 947
## [1] 948
## [1] 949
## [1] 950
## [1] 951
## [1] 952
## [1] 953
## [1] 954
## [1] 955
## [1] 956
## [1] 957
## [1] 958
## [1] 959
## [1] 960
## [1] 961
## [1] 962
## [1] 963
## [1] 964
## [1] 965
## [1] 966
## [1] 967
## [1] 968
## [1] 969
## [1] 970
## [1] 971
## [1] 972
## [1] 973
## [1] 974
## [1] 975
## [1] 976
## [1] 977
## [1] 978
## [1] 979
## [1] 980
## [1] 981
## [1] 982
## [1] 983
## [1] 984
## [1] 985
## [1] 986
## [1] 987
## [1] 988
## [1] 989
## [1] 990
## [1] 991
## [1] 992
## [1] 993
## [1] 994
## [1] 995
## [1] 996
## [1] 997
## [1] 998
## [1] 999
## [1] 1000

Examining Model Parameters

#####################################################
########### Examining Model Parameters ##############
#####################################################

par(mfrow=c(2,2))
hist(tausq.samp,main="tausq")
hist(mu0.samp,main="mu0")
hist(mu.samp[,1],main="mu group 1")
hist(mu.samp[,2],main="mu group 2")

par(mfrow=c(1,1))
boxplot(y)

Test posterior probabilities of gene expression for genes

# posterior probability group 5 has greater mean than group 6
postprob = sum(mu.samp[,5] > mu.samp[,6])/numsamp
postprob
## [1] 0.582
# posterior probability group 2 has greater mean than group 1
postprob = sum(mu.samp[,2] > mu.samp[,1])/numsamp
postprob
## [1] 0.019

Examining Shrinkage Graphically

#####################################################
######### Examining Shrinkage Graphically ###########
#####################################################

datameans = apply(y,2,mean)
postmeans = apply(mu.samp,2,mean)
mu0.mean = mean(mu0.samp)

par(mfrow=c(1,1))
plot(1:12,datameans,main="Shrinkage of Normal Means",pch=19)
abline(h=mu0.mean,col=4,lwd=2)
points(1:12,postmeans,pch=19,col=2)
legend(8,2,c("Data Mean","Post Mean","Mu0"),pch=19,col=c(1,2,4))

Posterior Predictive Sampling

#####################################################
######### Posterior Predictive Sampling #############
#####################################################

## sampling distribution of new observation 
## from a currently existing group

ystar.group1 = rep(NA,numsamp)
ystar.group2 = rep(NA,numsamp)
for (i in 1:numsamp){
    ystar.group1[i] = rnorm(1,mean=mu.samp[i,1],sd=sqrt(truesigsq))
    ystar.group2[i] = rnorm(1,mean=mu.samp[i,2],sd=sqrt(truesigsq))
}

par(mfrow=c(2,1))
xmin = min(c(ystar.group1,ystar.group2))
xmax = max(c(ystar.group1,ystar.group2))
hist(ystar.group1,main="Group 1 New Obs",xlim=c(xmin,xmax))
hist(ystar.group2,main="Group 2 New Obs",xlim=c(xmin,xmax))

## sampling distribution of new observation
## from an entirely new group

ystar.newgroup = rep(NA,numsamp)
for (i in 1:numsamp){
    mu.newgroup = rnorm(1,mean=mu0.samp[i],sd=sqrt(tausq.samp[i]))
    ystar.newgroup[i] = rnorm(1,mean=mu.newgroup,sd=sqrt(truesigsq))
}

par(mfrow=c(3,1))
xmin = min(c(ystar.group1,ystar.group2,ystar.newgroup))
xmax = max(c(ystar.group1,ystar.group2,ystar.newgroup))
hist(ystar.group1,main="Group 1 New Obs",xlim=c(xmin,xmax))
hist(ystar.group2,main="Group 2 New Obs",xlim=c(xmin,xmax))
hist(ystar.newgroup,main="New Group New Obs",xlim=c(xmin,xmax))