rasmusab / bayesian_first_aid Goto Github PK
View Code? Open in Web Editor NEWInside every classical test there is a Bayesian model trying to get out.
Inside every classical test there is a Bayesian model trying to get out.
Are you going to add things on the to-do list soon? bayes.var.test and bayes.wilcox.test would be pretty useful for some training I want to do.:-)
Hello,
I was trying to install the bayesian first aid package under R 3.1 running on OSX 10.7. I followed the instructions on the webpage for installing from github. Everything went well until the vignette. I let it run for about 2 hours before I killed the R session. Is there a way to install without the vignette?
Thanks!
OSX 10.7.5
R 3.1
BayesianFirstAid 0.1
JAGS 3.3.0 for OS X
I am playing (err serious work) around with some data for my PhD candidacy exams. I am running this command:
fit_meta <- bayes.cor.test( ~ KS.OTUs + Boron, data = secNDmeta)
Error in jags.model(textConnection(model_string), data = data, inits = inits, :
Error in node sigma[1]
Unobserved node inconsistent with observed parents at initialization
My data set looks like this (head version from R): Sorry the tab separated is screwy here:
AD.OTUs KS.OTUs X16s.OTUs ActiveCarbon Aluminum Antimonies Arsenic Barium Beryllium
AZ03 2263 2098 1070 168.43 11.01 0.47 0.06 20.40 0.47
AZ04 1739 1011 655 -6.67 20.36 0.20 0.08 14.71 0.12
AZ05 1060 733 911 7.63 7.49 0.00 0.31 23.28 0.15
AZ06 1992 1839 1027 172.00 10.92 0.34 0.05 22.94 0.29
AZ10 5129 3035 794 68.37 16.75 0.33 0.16 17.65 0.18
AZ11 1269 1427 687 -22.74 12.38 0.00 0.06 41.28 0.06
Boron Cadmium Calcium Chromium Cobalt Copper Iron KMNO4 LOI Lead Lithium Magnesium
AZ03 0.18 0.03 1687.48 0.00 0.02 0.29 0.73 0.51 2.91 0.28 0.11 298.06
AZ04 0.34 0.04 4545.20 0.00 0.00 1.08 0.26 0.53 0.95 0.09 0.14 157.88
AZ05 0.35 0.03 13526.16 0.02 0.03 0.33 0.80 0.53 2.50 0.08 0.19 148.07
AZ06 0.38 0.09 2718.54 0.00 0.01 0.47 0.43 0.51 2.53 0.36 0.07 163.69
AZ10 1.06 0.05 4283.93 0.01 0.01 0.87 0.50 0.52 1.25 0.36 0.35 510.39
AZ11 0.32 0.03 10311.61 0.00 0.01 0.26 0.73 0.53 0.78 0.24 0.09 155.72
Manganese Mercury Moisture Molybdenum Nickel Organic_Matter PMN Phosphorus Potassium
AZ03 4.32 0 1.09 0.00 0.06 1.81 10.77 7.07 223.27
AZ04 2.02 0 0.98 0.00 0.03 0.44 1.10 15.17 368.55
AZ05 8.92 0 0.40 0.00 0.05 1.52 3.81 10.28 146.51
AZ06 6.80 0 1.13 0.01 0.04 1.54 6.65 18.27 321.66
AZ10 7.50 0 1.64 0.00 0.08 0.65 4.61 29.12 415.41
AZ11 79.75 0 0.98 0.01 0.02 0.31 0.23 11.28 227.11
Selenium Silicon Sodium Strontium Sulfur Thallium Vanadium Zinc pH Altitude..m.
AZ03 0.02 29.76 438.70 10.39 4.05 0.04 0.01 0.72 6.73 1845.56
AZ04 0.00 41.78 500.13 12.51 8.19 0.03 0.01 0.25 8.31 1549.60
AZ05 0.00 17.51 1077.27 21.50 36.91 0.00 0.09 0.33 8.05 1536.19
AZ06 0.00 36.70 449.93 8.12 4.98 0.04 0.01 0.98 7.52 1472.49
AZ10 0.00 59.25 3232.91 25.01 70.45 0.02 0.02 0.26 8.07 1109.47
AZ11 0.00 18.60 610.42 0.00 11.33 0.00 0.00 0.33 8.41 1059.18
Latitude Longitude Geotype
AZ03 34.21755 -109.2290 Desert
AZ04 34.53717 -110.0940 Desert
AZ05 34.56431 -110.3824 Desert
AZ06 35.04198 -110.3869 Desert
AZ10 34.39214 -111.4503 Desert
AZ11 34.38891 -111.4580 Dry Forest
I did some quick google searching and found this error was common in older versions of JAGS but not in the new one. I ran another geochemical dataset using the same command early and had no problems.
I did find this:
As a side note, it is helpful in JAGS to provide initial values for the incompletely observed occupancy state z that are consistent with observed presences, as provided in this example with zinit. In other words if x(j,i,t,k)=1, provide an intial value of 1 for z(j,i,t). Unlike WinBUGS and OpenBUGS, if you do not do this, you’ll often (but not always) encounter an error message
Over at http://mbjoseph.github.io/blog/2013/02/24/com_occ/
But that doesn't seem to be quite what is wrong here.
Do you have any suggestions? Thank you for your time and help!
ara
On line 182 in bayes_prop_test.R, I believe the parameter cred_mass = cred.mass should be passed to mcmc_stats so that the specified HDI is computed for the theta differences.
Example current output:
bayes.prop.test(x = c(34, 54), n = c(400, 500), cred.mass = .9)$stats
mean sd HDI% comp HDIlo HDIup
theta[1] 0.08717148 0.01400703 90 0.5 0.06398359 0.10948452
theta[2] 0.10958137 0.01388894 90 0.5 0.08678420 0.13215677
x_pred[1] 34.89780000 7.88862708 90 0.5 22.00000000 47.00000000
x_pred[2] 54.80966667 9.98124649 90 0.5 38.00000000 70.00000000
theta_diff[1,2] -0.02240990 0.01966577 95 0.0 -0.06066122 0.01670726
Hi @rasmusab ,
I have weighted data that looks like this:
df <- data.frame(Year = 2012:2016,
Sub_A = c(221.0317, 310.5392, 247.7214, 184.7109, 243.2550),
Sub_B = c(276.4308, 322.6342, 227.5891, 184.3510, 188.2100),
Sub_C = c(127.01029, 99.89520, 70.19430, 73.02120, 76.79611),
Sub_D = c(443.7368, 411.5475, 346.2134, 390.3476, 375.8736),
Sub_E = c(236.0558, 258.3134, 199.3459, 164.9189, 222.2181),
Sub_F = c(615.1197, 591.3498, 482.6987, 450.9629, 400.5997),
Total = c(1253.000, 1339.751, 1033.000, 1051.000, 1148.000))
I'd like to compare the proportions for each subgroup from year to year, as per the following:
library(dplyr)
mutate(df, Sub_A_prop = Sub_A / Total,
Sub_B_prop = Sub_B / Total,
Sub_C_prop = Sub_C / Total,
Sub_D_prop = Sub_D / Total,
Sub_E_prop = Sub_E / Total,
Sub_F_prop = Sub_F / Total)
Due to the weighting effect on having decimals, I get the following error when I run bayes.prop.test
:
Failed check for discrete-valued parameters in distribution dbin
Is there a workaround? Also, is there the possibility of running multiple (3+) comparisons at the same time, or am I constrained to run dual comparisons one at time?
Thanks!
Phil
Hi,
I am trying to run a model but getting an error message. Here is my entire code. Could anyone please help me to understand why I am getting this error message. Thanks in advance
library(R2WinBUGS);
library(R2jags);
library(coda) ;
library(boot);
library(R2jags);
dat<-matrix(c(
6 ,5 , 1 ,1 , 1, 1,
5 ,1 , 1 ,3 , 1, 1,
6 ,5 , 6 ,1 , 1, 1,
6 ,6 , 3 ,1 , 1, 1,
6 ,6 , 6 ,1 , 2, 1,
6 ,6 , 6 ,1 , 1, 1,
6 ,6 , 6 ,6 , 6, 1,
6 ,6 , 6 ,1 , 1, 1,
6 ,6 , 6 ,6 , 6, 5,
6 ,6 , 6 ,1 , 1, 1,
6 ,6 , 1 ,2 , 1, 1,
6 ,6 , 6 ,6 , 4, 1), ncol=6)
N=nrow(dat) #index for rows
I=ncol(dat) #index for columns
rater=c(1,2,3,1,2,3)
item=c(1,1,1,2,2,2)
K <- apply(dat, 2, max)
rater1 <- function() {
for (i in 1:N) {
for (j in 1:I) {
dat[i,j] ~ dcat(prob[i,j,1:K[j]])
}
theta[i] ~ dnorm(0,1)
for(j in 1:I) {
for(k in 1:(K[j]-1)) {
logit( P[i,j,k] ) <- theta[i]-B[item[j]] -C[item[j],k]-D[rater[j]]
}
P[i, j,K[j]] <- 1
}
for(j in 1:I) {
prob[i, j, 1] <- P[i, j, 1]
for (k in 2:(K[j])){
prob[i, j, k] <- P[i, j, k] - P[i, j, k-1]
}
}
}
for(j in 1:2) {
B[j] ~ dnorm(0,4)
}
for (j in 1:3) {
D[j] ~ dnorm(0,4)
}
for (j in 1:2) {
for(k in 1:(K[j])) {
C[j,k] ~ dnorm(0,1)
}
}
}
write.model(rater1, con = "rater.bug", digits = 5)
data=list("dat", "N", "I","rater","item", "K")
#R2jugs
rjags=jags(data=data, inits = NULL, n.chains = 1, n.iter = 1500,
n.thin = 1, model.file = "rater.bug", n.burnin = 300,DIC = T,
parameters.to.save=c("B","C","theta","D"))
Heading says "Bayesian First Aid propotion test" when it probably should be "Bayesian First Aid proportion test"
n_right_leaners <- c(4, 8)
n_respondents <- c(16, 16)
bayes.prop.test(n_right_leaners, n_respondents, n.iter = 15000)
Bayesian First Aid **propotion** test
data: n_right_leaners out of n_respondents
number of successes: 4, 8
number of trials: 16, 16
Estimated relative frequency of success [95% credible interval]:
Group 1: 0.27 [0.093, 0.48]
Group 2: 0.50 [0.28, 0.72]
Estimated group difference (Group 1 - Group 2):
-0.23 [-0.52, 0.074]
The relative frequency of success is larger for Group 1 by a probability
of 0.079 and larger for Group 2 by a probability of 0.921 .
I am new to R, so apologies if this has a simple answer. Let's say I've fitted a model... rather than printing out the summary information using summary(fit)
, is it possible to assign the summary stats (i.e. mean sd HDIlo HDIup etc) to a data frame? I'm attempting to plot posterior summary information in RMarkdown/Knitr, so just wondered how I can extract the info into a variable.
PS, great repo
Thanks!
I do not know the source of my problem. Does the following mean anything to you? Is so can you suggest a troubleshooting avenue or two that I can go down.
install_github("rasmusab/bayesian_first_aid")
has always worked in the past but now I get this
> install_github("rasmusab/bayesian_first_aid")
Downloading GitHub repo rasmusab/bayesian_first_aid@master
from URL https://api.github.com/repos/rasmusab/bayesian_first_aid/zipball/master
Installing BayesianFirstAid
"C:/PROGRA1/R/R-331.0/bin/i386/R" --no-site-file --no-environ --no-save
--no-restore --quiet CMD INSTALL
"C:/Users/Farrel/AppData/Local/Temp/RtmpWEoqvP/devtools2e9c79894018/rasmusab-bayesian_first_aid-6c7f844"
--library="G:/Farrel Data/Documents/R/win-library/3.3" --install-tests
More of a suggestion than an issue. I had difficulty finding documentation on making rjags/JAGs play nicely. It might probably worth adding a link to this gist https://gist.github.com/casallas/8411082 in the README.
love the lib btw!
Do you have any plans for venturing into linear regression, multiple linear regression, logistic regression and multiple logistic regression and survival analysis. I personally do not know how to code it but I want to encourage you or anyone else you can coopt to write with you.
I am unable to install the package on R 3.6.3. I get the following error. I am using Google Colab's R notebook (colab.fan/r)
devtools::install_github("rasmusab/bayesian_first_aid")
Downloading GitHub repo rasmusab/bayesian_first_aid@master
rjags (NA -> 4-10 ) [CRAN]
coda (NA -> 0.19-3) [CRAN]
Skipping 1 packages not available: mnormt
Installing 3 packages: rjags, coda, mnormt
Installing packages into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
Error: Failed to install 'BayesianFirstAid' from GitHub:
(converted from warning) package ‘mnormt’ is not available (for R version 3.6.3)
Traceback:
In RStudio at least, calling plot()
or diagnostics()
on a fitted model causes any subsequent base graphics to be drawn inside one of the subplot regions used by BayesianFirstAid
.
Example:
b.t = bayes.t.test(subject.means.conf.foil$md)
plot(b.t)
hist(subject.means.conf.foil$md)
The problem is easily solved by calling par(mfrow=c(1,1))
in between plots, as far as I can tell, but this could maybe be incorporated into the package?
Hi,
I was wondering whether it is possible to set a seed that generates reproducible results? I noticed that when running tests, I get slightly different results each time.
My understanding is that the package calls rjags. So the seed would have to be set at that step. Based on the rjags manual, it should be possible to set a RNG (Random number generators) for the inits
parameter in the jags.model
function.
Thanks,
Respectfully request Windows binary version some day ...
This just started happening:
`> x <- c(122, 99)
n <- c(4923, 5001)
bayes.prop.test(x, n)
Error in get(as.character(FUN), mode = "function", envir = envir) :
object 'home_state' of mode 'function' was not found`
I have the following packages loaded:
`> sessionInfo()
R version 3.1.1 Patched (2014-08-27 r66482)
Platform: x86_64-apple-darwin10.8.0 (64-bit)
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] BayesianFirstAid_0.1 rjags_3-13 coda_0.16-1
[4] lattice_0.20-29 lazyeval_0.1.9 stringr_0.6.2
[7] reshape2_1.4 PerformanceAnalytics_1.4.3541 magrittr_1.0.1
[10] wf_1.0 timeDate_3010.98 quantmod_0.4-0
[13] TTR_0.22-0 Defaults_1.1-1 xts_0.9-7
[16] zoo_1.7-11 rjson_0.2.14 testthat_0.8.1
[19] roxygen2_4.0.2 devtools_1.5 RCurl_1.95-4.3
[22] bitops_1.0-6 lubridate_1.3.3 assertthat_0.1
[25] scales_0.2.4 ggthemes_1.7.0 tidyr_0.1
[28] ggplot2_1.0.0 dplyr_0.3.0.2 RPostgreSQL_0.4
[31] DBI_0.3.0
loaded via a namespace (and not attached):
[1] colorspace_1.2-4 digest_0.6.4 evaluate_0.5.5 grid_3.1.1 gtable_0.1.2
[6] httr_0.5 labeling_0.3 MASS_7.3-34 memoise_0.2.1 mnormt_1.5-1
[11] munsell_0.4.2 parallel_3.1.1 plyr_1.8.1 proto_0.3-10 RColorBrewer_1.0-5
[16] Rcpp_0.11.2 tools_3.1.1 whisker_0.3-2 `
This does not happen with a fresh instance of R, so probably a masking problem.
I was using BFA for a quick analysis of a Candy Japan A/B test, replacing prop.test
with a Bayesian test, and then demonstrating the difference that an informative prior on small effect sizes (in this case dbeta(900, 1407)
) makes. The model.code
functionality is very nice and convenient, but it still requires something like 18 lines to extract and copy-paste and isn't very friendly to the reader. It would be nice if I could do something like
fit1 <- bayes.prop.test(c(175, 168), c(442, 439))
# with informative prior:
fit2 <- bayes.prop.test(c(175, 168), c(442, 439), prior="dbeta(900,1407)")
And a prior
parameter gets substituted in the JAGS model code. (I peeked at the implementations but it looks like it bottoms out at something called jags_binom_test
which isn't available.)
By the way, there seems to be a typo in bayes.prop.test
where it prints "Bayesian First Aid poportion test" instead of "Bayesian First Aid proportion test".
Hi,
I'm trying to install bayesfirstaid as instructed with rtools and JAGS installed and this is what I get:
devtools::install_github("rasmusab/bayesian_first_aid")
Downloading GitHub repo rasmusab/bayesian_first_aid@HEAD
Skipping 1 packages not available: mnormt
√ checking for file 'C:\Users\dk3434\AppData\Local\Temp\RtmpA5fVvn\remotes51b4555f4e43\rasmusab-bayesian_first_aid-d80c0fd/DESCRIPTION' (507ms)
Hi,
I've just discovered your nice R package because I was searching for a Bayesian estimation applicable for data for which you would normally use a paired t.test.
My data consists of acceptability ratings given my the same subjects in two differing experimental conditions. (The acceptability ratings might actually be better treated as ordinal data but I'm treating them as metrical).
Using a default
t.test(d$condition1, d$condition2, paired=TRUE)
gives this result:
Paired t-test
data: d$condition1 and d$condition2
t = -9.3649, df = 407, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.7561958 -0.4938042
sample estimates:
mean of the differences
-0.625
As I read this (and the data also looks like this), this is quite good evidence that people rated stimuli in condition1 lower than in condition2.
However, if I'm using
bayes.t.test(d$condition1, d$condition2, paired= TRUE)
Bayesian estimation supersedes the t test (BEST) - paired samples
data: d$condition1 and d$condition2, n = 408
Estimates [95% credible interval]
mean paired difference: -5.6e-07 [-0.00012, 0.00013]
sd of the paired differences: 0.0014 [0.0013, 0.0014]
The mean difference is more than 0 by a probability of 0.497
and less than 0 by a probability of 0.503
These results tell a complete opposite: Apparently, there is no difference under the two conditions. What looks strange to me is that the estimated mean and SD do not match the mean and SD of the data at all (which is M = 0.625 and SD = 1.35)
Since I'm just discovering Bayesian statistics, I wonder if I'm something wrong or if there is some bug in your R package.
You may find the data for the next 2 weeks under this link.
Hi,
First, thanks so much for producing this library, it's great.
I'm using the code and producing an html file via Knitr in RStudio. When I activate the library, using:
library(BayesianFirstAid, warn.conflicts=FALSE, verbose=FALSE, quietly = TRUE)
I get the following output in my html file:
I'm not very familiar with html, I'm sure I could edit that out, but it would be nice to just get it to load silently.
Thanks for any help you can provide.
Steven
I do not know the source of my problem. Does the following mean anything to you? If so can you suggest a troubleshooting avenue or two that I can go down.
install_github("rasmusab/bayesian_first_aid")
has always worked in the past but now I get this
> install_github("rasmusab/bayesian_first_aid")
Downloading GitHub repo rasmusab/bayesian_first_aid@master from URL https://api.github.com/repos/rasmusab/bayesian_first_aid/zipball/master Installing BayesianFirstAid "C:/PROGRA~1/R/R-33~1.0/bin/i386/R" --no-site-file --no-environ --no-save \ --no-restore --quiet CMD INSTALL \ "C:/Users/Farrel/AppData/Local/Temp/RtmpWEoqvP/devtools2e9c79894018/rasmusab-bayesian_first_aid-6c7f844" \ --library="G:/Farrel Data/Documents/R/win-library/3.3" --install-tests
* installing *source* package 'BayesianFirstAid' ... ** R Warning: unable to re-encode 'bayes_prop_test.R' line 8 ** inst ** preparing package for lazy loading Error: Command failed (5)
Hi,
I tried to install the package but got this issue.
devtools::install_github("rasmusab/bayesian_first_aid")
Error: Failed to install 'unknown package' from GitHub:
Line starting 'Roxyg ...' is malformed!
remotes::install_github("rasmusab/bayesian_first_aid")
Error: Failed to install 'unknown package' from GitHub:
Line starting 'Roxyg ...' is malformed!
devtools::install_github("rasmusab/bayesian_first_aid", upgrade = 'never')
Error: Failed to install 'unknown package' from GitHub:
Line starting 'Roxyg ...' is malformed!
Anyone knows how to solve it? Thank you!
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.