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"
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"
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
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)
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))