library(rstan)
library(brms)
library(yaml)
library(bayesplot)
library(gtsummary)
library(dplyr)
library(blastula)
## tidyberse
library(flagr)
library(sjPlot)
library(insight)
library(httr)
library(ggthemes)
library(tidybayes)
library(ggplot2)
library(gt)
library(scales)
library(cowplot)
library(patchwork)
library(ghibli)
library(broom)
library(stringr)
library(StanHeaders)
library(tidyr)
#addline_formatの関数：これでラベルを2段組みにできる
addline_format <- function(x,...){
gsub('\\s','\n',x)
}
#解析中のメモリの負担を減らす
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#++++++++++++++++++++++++++分析+++++++++++++++++++++++++++++++++++++++++++++++++++
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
# データの読み込み
dd <- read.csv("tt.csv")
tt0 <- brm(
formula = cellcount ~ condition ,
family = negbinomial(),#poissonよりこっちやな
data = dd,
seed = 1,
iter = 9000,
warmup = 6000,
chains = 4,
thin = 1,
control = list(adapt_delta = 0.99, max_treedepth = 15))
## 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.2-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-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.2-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.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2-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.2-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.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2-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
plot(tt0)
mcmc_plot(tt0,type = "rhat")
summary(tt0)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: cellcount ~ condition
## Data: dd (Number of observations: 32)
## Draws: 4 chains, each with iter = 9000; warmup = 6000; thin = 1;
## total post-warmup draws = 12000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 4.72 0.18 4.38 5.10 1.00 8393 6553
## conditionlive 0.37 0.25 -0.14 0.87 1.00 8296 6565
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 2.11 0.52 1.22 3.26 1.00 8040 6525
##
## 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).
brms::pp_check(tt0)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#++++++++++++++++分析データを作図のために調節+++++++++++++++++++++++++
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ndd <- crossing(nesting(dd$condition,dd$nodo))
ndd #変数の組み合わせを確認するだけ
## # A tibble: 8 × 2
## `dd$condition` `dd$nodo`
## <chr> <chr>
## 1 death n0
## 2 death n1
## 3 death n2.5
## 4 death n5
## 5 live n0
## 6 live n1
## 7 live n2.5
## 8 live n5
# A tibble: 8 × 2
# `dd$condition` `dd$nodo`
# <chr> <chr>
#1 death n0
#2 death n1
#3 death n2.5
#4 death n5
#5 live n0
#6 live n1
#7 live n2.5
#8 live n5
df <- data.frame(matrix(vector(),8,5,dimnames=list(c(), c("condition","nodo","Mean","low","high"))),stringsAsFactors=F)
df$nodo <- c("n0","n1","n2.5","n5","n0","n1","n2.5","n5")
df$condition <- c("death","death","death","death","live","live","live","live")
pisize <- c(0.025, 0.975)
######################################
mod <- tt0 # cellcount ~ condition
######################################
df[1,3:5] <- fitted(mod, newdata = tibble(nodo= 'n0', condition="death") ,allow_new_levels = TRUE,probs = pisize) [,c(1,3,4)]
df[2,3:5] <- fitted(mod, newdata = tibble(nodo= 'n1',condition="death") ,allow_new_levels = TRUE,probs = pisize) [,c(1,3,4)]
df[3,3:5] <- fitted(mod, newdata = tibble(nodo= 'n2.5',condition="death") ,allow_new_levels = TRUE,probs = pisize) [,c(1,3,4)]
df[4,3:5] <- fitted(mod, newdata = tibble(nodo= 'n5',condition="death") ,allow_new_levels = TRUE,probs = pisize) [,c(1,3,4)]
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
df[5,3:5] <- fitted(mod, newdata = tibble(nodo= 'n0',condition="live") ,allow_new_levels = TRUE,probs = pisize) [,c(1,3,4)]
df[6,3:5] <- fitted(mod, newdata = tibble(nodo= 'n1',condition="live") ,allow_new_levels = TRUE,probs = pisize) [,c(1,3,4)]
df[7,3:5] <- fitted(mod, newdata = tibble(nodo= 'n2.5',condition="live"),allow_new_levels =TRUE,probs = pisize) [,c(1,3,4)]
df[8,3:5] <- fitted(mod, newdata = tibble(nodo= 'n5',condition="live"),allow_new_levels = TRUE,probs = pisize) [,c(1,3,4)]
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#++++++++++++++++++++++++++分析データの作図+++++++++++++++++++++++++++++++++++++++
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
m<-"royalblue" ; f<-"grey40"
pp <- ggplot(data=df,aes(x=as.factor(nodo),y=Mean, ymin=low,ymax=high,
color=as.factor(condition),shape= as.factor(condition))) +
# geom_line(data=df,aes(y = Mean, group = as.factor(condition)),
# linewidth = 1,
# color = "grey70",
# linetype = "dotted") +
geom_pointrange(data=df,position=position_dodge(width=0.3),fatten = 5,size = 1) +
scale_shape_manual(values = c(16, 17),# 16=◍、17=▲、1=◯、2＝△
name="",
labels=c("death","live"))+
# scale_y_continuous(breaks=seq(12,350,12))+
scale_color_manual(values=c(m,f),
name ="",
labels=c("death","live")) +
ggtitle("") + xlab("nodo") + ylab("countcell")+ #
#累積再配偶回数と出産間隔の関係
#facet_wrap(~Res,strip.position = "bottom")+
theme_bw() +
theme(strip.placement = "outside",strip.background = element_blank(),plot.title = element_text(hjust = 0.5),panel.grid = element_blank(),text = element_text(size =16,face = ))+
theme(legend.position = "bottom")+
theme(axis.title.x = element_text(size =18,face="bold"))+
theme(axis.title.y = element_text(size =18,face="bold"))+
theme(axis.text.y = element_text(size =15))+
theme(axis.text.x = element_text(size =18,face="bold"))+
theme(strip.text.x = element_text(size =17))+
theme(legend.text = element_text(size =17))+
theme(legend.title = element_text(size =17))+
theme(legend.box = "vertical")
pp # 一応完成
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#++++++++++++++++++++++++++++生データ+++++++++++++++++++++++++++++++++++++++++++++
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#色指定はここですると後々楽
m<-"royalblue" ; f<-"grey40"
s <-"black"
tt <- ggplot(dd, aes(x = nodo, y = cellcount, fill = condition)) +
geom_bar(position = "dodge",stat = "identity",alpha=0.7,width = 0.5)+
# もしx軸の数値を調節したかったら、このコード使って（左から希望の最低値、最高値、〜ごとの区切り、この例だと２４ごとにカウントされる） scale_x_continuous(breaks=seq(12,350,24))+
scale_fill_manual(values=c(m,f),
name=addline_format("cell status"),
labels=c("live","death"))+
ggtitle("cell count results by concentration")+
xlab("nodo(add a unit)") +
ylab("cell count(add a unit)") +
#文字サイズとか細かい設定（任意）
theme(axis.title.x = element_text(size =15))+
theme(axis.title.y = element_text(size =15))+
theme(axis.text.y = element_text(size =15))+
theme(axis.text.x = element_text(size =15))+
theme(strip.text.x = element_text(size =15))+
theme(legend.text = element_text(size =8))+
theme_bw() +
theme(legend.position="right")+
theme(strip.placement = "outside",strip.background = element_blank(),plot.title = element_text(hjust = 0.5),panel.grid = element_blank(),text = element_text(size =16,face = ))
# facet_wrap(~celltype,ncol = 3)
tt
