/
summary_stat.R
158 lines (144 loc) · 7.74 KB
/
summary_stat.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#' Summary statistics from posterior of MixSIAR model
#'
#' \code{summary_stat} prints and saves summary statistics
#'
#' @param combined list, output from \code{\link{combine_sources}} function
#' @param toprint vector, which parameters to print? Options are: \code{"p"} to print stats for proportions only (all factors),
#' \code{"global"} to only print overall proportions, \code{"fac1"} to only print factor 1 proportions, \code{"fac2"} to print
#' factor 2 proportions. Set = \code{"epsilon"} to print multiplicative error term only.
#' Default = \code{"all"}, prints stats for all model parameters.
#' @param groupby character, group stats by "factor" or "source"? I.e. in wolves example, group proportions by Region 1, Region 2, Region 3
#' (\code{groupby="factor"}) vs. Deer, Marine Mammals, Salmon (\code{groupby="source"}). Currently only \code{"factor"} is implemented.
#' @param meanSD \code{TRUE/FALSE}, print mean and SD for the parameters?
#' @param quantiles vector, which quantiles to print. Default = \code{c(0.025, 0.25, 0.5, 0.75, 0.975)}.
#' @param savetxt \code{TRUE/FALSE}, save results as .txt file (in working directory)?
#' @param filename character, file name to save results as (\code{.txt} will be appended automatically)
#'
#' @return NULL
#' @export
#' @seealso \code{\link{combine_sources}} and \code{\link{plot_intervals}}
#' @examples
#' \dontrun{
#' # first run mantis shrimp example
#' # combine 6 sources into 2 groups of interest (hard-shelled vs. soft-bodied)
#' # 'hard' = 'clam' + 'crab' + 'snail' # group 1 = hard-shelled prey
#' # 'soft' = 'alphworm' + 'brittlestar' + 'fish' # group 2 = soft-bodied prey
#' combined <- combine_sources(jags.1, mix, source, alpha.prior=alpha,
#' groups=list(hard=c("clam","crab","snail"), soft=c("alphworm","brittlestar","fish")))
#'
#' summary_stat(combined)
#' summary_stat(combined, savetxt=FALSE)
#' summary_stat(combined, meanSD=FALSE)
#' summary_stat(combined, quantiles=c(.05,.5,.95))
#' summary_stat(combined, toprint="fac1")
#' summary_stat(combined, toprint="p")
#' summary_stat(combined, toprint="global")
#' }
summary_stat <- function(combined, toprint="all", groupby="factor", meanSD=TRUE, quantiles=c(.025,.25,.5,.75,.975), savetxt = TRUE, filename="summary_statistics"){
# Error check for groupby
if(!(groupby %in% c("factor","source"))){
stop(paste0("groupby must be 'factor' or 'source'.
See ?summary_stat"))
}
if(groupby=="source") warning("My apologies... groupby='source' not implemented yet.")
# Error check for toprint
if(!(toprint %in% c("all","p","global","fac1","fac2","epsilon"))){
stop(paste0("toprint options are 'all', 'p', 'global', 'fac1', 'fac2', 'epsilon'.
See ?summary_stat"))
}
if(combined$mix$n.fe != 0 & toprint=="global"){
stop(paste0("toprint='global' but no global proportions in model.
You have fixed effects in your model, so intercept terms are meaningless.
Set toprint='fac1' to print factor 1 proportions,
or toprint='fac2' to print factor 2 proportions."))
}
post <- combined$post
n.sources <- combined$source.new$n.sources
source_names <- combined$source.new$source_names
# easier - just get the proportions from post as new objects
p.global <- post[,grep("p.global",colnames(post))]
p.fac1 <- post[,grep("p.fac1",colnames(post))]
p.fac2 <- post[,grep("p.fac2",colnames(post))]
p.both <- post[,grep("p.both",colnames(post))]
# groupby="source" -- need to reorder p.fac1, p.fac2, p.both
stats <- NULL; eps_stats <- NULL; eps_labels <- NULL;
sig_labels <- NULL; fac1_labels <- NULL; fac2_labels <- NULL; sig_stats <- NULL;
getQuant <- function(x) quantile(x, probs=quantiles)
getMeanSD <- function(x) cbind(round(apply(x,2,mean),3),round(apply(x,2,sd),3))
if(combined$mix$n.fe == 0 & (toprint %in% c("all","p","global"))){
global_quants <- t(round(apply(p.global,2,getQuant),3))
if(meanSD){ global_means <- getMeanSD(p.global)} else {global_means <- NULL}
stats <- cbind(global_means, global_quants)
global_labels <- paste("p.",source_names,sep="")
rownames(stats) <- global_labels
}
if(combined$mix$n.effects > 0 & combined$mix$n.fe != 2 & (toprint %in% c("all","p","fac1"))){
if(groupby=="factor"){
fac1_quants <- t(round(apply(p.fac1,2,getQuant),3))
if(meanSD){ fac1_means <- getMeanSD(p.fac1)} else {fac1_means <- NULL}
fac1_stats <- cbind(fac1_means,fac1_quants)
fac1_labels <- paste0("p.",do.call(paste, c(expand.grid(source_names,combined$mix$FAC[[1]]$labels),sep=".")))
rownames(fac1_stats) <- fac1_labels
stats <- rbind(stats,fac1_stats)
}
if(combined$mix$FAC[[1]]$re & (toprint %in% c("all","fac1"))){
sig_stats <- cbind(getMeanSD(combined$jags.1$BUGSoutput$sims.list$fac1.sig),t(round(apply(combined$jags.1$BUGSoutput$sims.list$fac1.sig,2,getQuant),3)))
sig_labels <- paste(combined$mix$FAC[[1]]$name,".SD",sep="")
}
}
if(combined$mix$n.re==2 & (toprint %in% c("all","p","fac2"))){
if(groupby=="factor"){
fac2_quants <- t(round(apply(p.fac2,2,getQuant),3))
if(meanSD){ fac2_means <- getMeanSD(p.fac2)} else {fac2_means <- NULL}
fac2_stats <- cbind(fac2_means,fac2_quants)
fac2_labels <- paste0("p.",do.call(paste, c(expand.grid(source_names,combined$mix$FAC[[2]]$labels),sep=".")))
rownames(fac2_stats) <- fac2_labels
stats <- rbind(stats,fac2_stats)
}
if(combined$mix$FAC[[2]]$re & (toprint %in% c("all","fac2"))){
sig_stats <- rbind(sig_stats,cbind(getMeanSD(combined$jags.1$BUGSoutput$sims.list$fac2.sig),t(round(apply(combined$jags.1$BUGSoutput$sims.list$fac2.sig,2,getQuant),3))))
sig_labels <- c(sig_labels,paste(combined$mix$FAC[[2]]$name,".SD",sep=""))
}
}
if(combined$mix$fere & (toprint %in% c("all","p","fac1","fac2"))){
fac2_quants <- t(round(apply(p.both,2,getQuant),3))
if(meanSD){ fac2_means <- getMeanSD(p.both)} else {fac2_means <- NULL}
fac2_labels <- paste0("p.",do.call(paste, c(expand.grid(source_names,combined$mix$FAC[[2]]$labels,combined$mix$FAC[[1]]$labels),sep=".")))
fac2_stats <- round(cbind(fac2_means,fac2_quants),3)
rownames(fac2_stats) <- fac2_labels
stats <- rbind(stats,fac2_stats)
if(combined$mix$FAC[[2]]$re & (toprint %in% c("all","fac2","fac1"))){
sig_stats <- rbind(sig_stats,cbind(getMeanSD(combined$jags.1$BUGSoutput$sims.list$fac2.sig),t(round(apply(combined$jags.1$BUGSoutput$sims.list$fac2.sig,2,getQuant),3))))
sig_labels <- c(sig_labels,paste(combined$mix$FAC[[2]]$name,".SD",sep=""))
}
}
# Add SD stats to the top of the summary
rownames(sig_stats) <- sig_labels
stats <- rbind(sig_stats,stats)
# Add epsilon (multiplicative error term) to stat summary
# Also plot posterior density
epsTF <- "resid.prop" %in% names(combined$jags.1$BUGSoutput$sims.list)
if(epsTF & (toprint %in% c("all","epsilon"))){
if(meanSD) eps_stats <- cbind(getMeanSD(combined$jags.1$BUGSoutput$sims.list$resid.prop),t(round(apply(combined$jags.1$BUGSoutput$sims.list$resid.prop,2,getQuant),3)))
if(!meanSD) eps_stats <- t(round(apply(combined$jags.1$BUGSoutput$sims.list$resid.prop,2,getQuant),3))
eps_labels <- paste0("Epsilon.", 1:combined$mix$n.iso)
rownames(eps_stats) <- eps_labels
stats <- rbind(eps_stats,stats)
}
if(meanSD) colnames(stats) <- c("Mean","SD", paste(round(100*quantiles, 1), "%", sep=""))
if(!meanSD) colnames(stats) <- paste(round(100*quantiles, 1), "%", sep="")
out1 <- capture.output(stats)
cat("
################################################################
# Summary Statistics
################################################################
",out1,sep="\n")
if(savetxt){ # svalue(summary_save)
mypath <- file.path(paste(getwd(),"/",filename,".txt",sep=""))
cat("
#################################################################
# Summary Statistics
#################################################################
",out1,sep="\n", file=mypath, append=TRUE)
}
}