Title: | Analysis of Scientific Evidence Using Bayesian and Likelihood Methods |
---|---|
Description: | Bayesian (and some likelihoodist) functions as alternatives to hypothesis-testing functions in R base using a user interface patterned after those of R's hypothesis testing functions. See McElreath (2016, ISBN: 978-1-4822-5344-3), Gelman and Hill (2007, ISBN: 0-521-68689-X) (new edition in preparation) and Albert (2009, ISBN: 978-0-387-71384-7) for good introductions to Bayesian analysis and Pawitan (2002, ISBN: 0-19-850765-8) for the Likelihood approach. The functions in the package also make extensive use of graphical displays for data exploration and model comparison. |
Authors: | Robert van Hulst |
Maintainer: | Robert van Hulst <[email protected]> |
License: | GPL (>= 2) |
Version: | 0.8.10 |
Built: | 2024-11-02 04:28:31 UTC |
Source: | https://github.com/cran/evidence |
The functions in this package include Bayesian and likelihood alternatives to the standard statistical hypothesis tests that form part of base R. Their aim is to provide a wider perspective on how statistical evidence can be analyzed than the usual hypothesis-testing one. In view of the increasing importance in science of Bayesian and likelihood inference a wider exposure to these alternatives has become overdue.
This package makes Bayesian and likelihood analyses of simple statistical problems as convenient as traditional frequentist ones are in R. In addition, it makes effective use of R's excellent plotting capabilities, and facilitates exploratory data analysis and an interactive approach to modeling. Both data exploration and model exploration are crucial in data analysis, and these are facilitated by an interactive and graphics-centered approach.
Package: | evidence |
Type: | Package |
Version: | 0.8.10 |
Date: | 2018-04-15 |
License: | GPL |
Robert van Hulst
Maintainer: <rvhulst.ubishops.ca>
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
Made-up data with easy numbers for practicing one-way anova by hand to understand how an anova works.
data(AOV1)
data(AOV1)
A data frame with 15 observations on the following 2 variables.
y
response
i
predictor, a factor with 3 levels
Note that the design is balanced.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(AOV1) summary(aov(y ~ i, data=AOV1))
data(AOV1) summary(aov(y ~ i, data=AOV1))
Made-up data with easy numbers for practicing one-way anova by hand to understand how an anova works.
data(AOV2)
data(AOV2)
A data frame with 22 observations on the following 2 variables:
y
response
i
predictor: a factor with 4 levels
Note that the design is unbalanced.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(AOV2, package) summary(aov(y ~ i, data=AOV2))
data(AOV2, package) summary(aov(y ~ i, data=AOV2))
The Physicians health study data cross-classified according to Infarct (heart attack or not) and Group (Placebo or Aspirin).
data(Aspirin)
data(Aspirin)
A 2 by 2 matrix of counts with row names:
Infarct:Yes and Infarct:No,
and column names:
Group:placebo and Group:aspirin.
Steering Committee of the Physicians' Health Study Research Group. 1989. Final report of the aspirin component of the ongoing Physicians' Health Study. N Engl J Med, 321:129–135.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
This function performs a standard Bayesian analysis of a single sample of a population presumably following a Normal distribution. Imprecise priors for the mean and the standard deviation are used.
B1Nmean(x, plotit = TRUE, hists = FALSE, pdf = FALSE)
B1Nmean(x, plotit = TRUE, hists = FALSE, pdf = FALSE)
x |
a vector of sample values |
plotit |
should the function produce plots? Defaults to TRUE. |
hists |
should histograms of the posterior distribution for the data with twenty posterior predictive histograms also be plotted? Defaults to FALSE. |
pdf |
should the histograms be saved as a pdf-file? Defaults to FALSE. |
none produced: text and graphical output are produced
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
## Not run: data(Fat) B1Nmean(Fat$Height) ## End(Not run)
## Not run: data(Fat) B1Nmean(Fat$Height) ## End(Not run)
This function performs a standard Bayesian analysis of a single sample of a population assumed to follow a Normal distribution. A Standard Improper Reference prior is assumed.
B1Nsir(x, r = 10000, alpha = 0.05)
B1Nsir(x, r = 10000, alpha = 0.05)
x |
a vector of sample values |
r |
the number of samples to be taken from the posterior distribution (defaults to 10000) |
alpha |
1 - level of credibility, so that for alpha = 0.05 (the default) credible intervals will have 95% credibility |
none returned; the function produces a plot of the posterior distribution and prints some statistics.
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(darwin) B1Nsir(darwin$difference)
data(darwin) B1Nsir(darwin$difference)
This function computes the posterior distribution of the binomial
probability when given the number of “successes” and the sample
size, as well as one of a choice of priors. A plot of the posterior
distribution is produced with the 95% credible interval of
.
B1prop(s, n, p = 0.5, alpha = 0.05, prior = c("uniform", "near_0.5", "not_near_0.5", "near_0", "near_1", "custom"), params = NULL)
B1prop(s, n, p = 0.5, alpha = 0.05, prior = c("uniform", "near_0.5", "not_near_0.5", "near_0", "near_1", "custom"), params = NULL)
s |
the number of sampling units with the feature |
n |
the number of sampling units examined |
p |
an optional hypothesized probability |
alpha |
1 - alpha is the desired level of credibility of a credible interval |
prior |
one of: "uniform", "near_0.5", "not_near_0.5", "near_0", "near_1", "custom", which are all beta distributions with appropriate parameter values. Note that if prior="custom" the following argument has to be supplied: |
params |
a vector with the a and b parameters of the custom beta prior |
the posterior probability
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
B1prop(13, 100, .1, prior="near_0")
B1prop(13, 100, .1, prior="near_0")
.
Provides a simple demonstration of how the posterior distribution improves as increasing amounts of data become available. A Binomial variable with a known parametric probability is sampled, and as increasing numbers of samples become available the posterior distribution is re-evaluated and plotted.
B1propSim(p, N = 100, prior = c("uniform", "near_0.5", "not_near_0.5", "near_0", "near_1"))
B1propSim(p, N = 100, prior = c("uniform", "near_0.5", "not_near_0.5", "near_0", "near_1"))
p |
the “real” binomial probability; if a number samller than 0 or one lager than 1 isentered the function will choose an arbitrary probability |
N |
the number of observations to accumulate |
prior |
one of: "uniform", "near_0.5", "not_near_0.5", "near_0", or "near_1". |
none returned; the function is run for the plot it produces.
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
B1propSim(p = 0.44, prior = "near_0.5")
B1propSim(p = 0.44, prior = "near_0.5")
Produces exploratory plots (boxplots and, if the sample sizes are equal), a quantile-quantile plot of the two samples. Also produces Bayesian posterior densities of the two sample means and of the difference between the means. The priors used are standard improper reference priors.
B2Nsir(formula, data, var.equal = TRUE, alpha = 0.05, plotit = TRUE, r = 10000)
B2Nsir(formula, data, var.equal = TRUE, alpha = 0.05, plotit = TRUE, r = 10000)
formula |
the standard formula interface: response ~ factor |
data |
a data.frame containing the response and the two-level factor |
var.equal |
if TRUE the group variances are assumed to be equal, if FALSE two separate group variances are estimated |
alpha |
1 - level of credibility, so that for alpha = 0.05 (the default) credible intervals will have 95% credibility |
plotit |
should plots be produced? |
r |
the number of samples from the posterior distribution; can usually be left at its default value of 10000 |
Note that in the first plot the second sub-plot is NOT a normality plot but a quantile-quantile plot that compares the observations in the two groups.
none returned; the function produces several plots and prints some statistics.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(bodytemp) B2Nsir(temperature ~ gender, bodytemp)
data(bodytemp) B2Nsir(temperature ~ gender, bodytemp)
This function computes the posterior distributions of the binomial
parameters and
when given the numbers of
“successes” and the sample sizes for the two samples. It uses uniform
priors. A plot of the posterior distributions of the two
's is
produced, and a plot of the posterior distribution of
with its 95% credible interval.
B2props(s, n, alpha = 0.05)
B2props(s, n, alpha = 0.05)
s |
a vector containing the 2 numbers of sampling units with the feature ("success") |
n |
a vector containing the 2 numbers of sampling units examined |
alpha |
1 - level of credibility, so that for alpha = 0.05 (the default) credible intervals will have 95% credibility |
None, the inferred difference between the probabilities and its 95% credible interval is calculated and several plots are produced
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
B2props(c(13, 22), c(78, 92))
B2props(c(13, 22), c(78, 92))
A 2 x 2 contingency table (in matrix form) is analyzed in a Bayesian way using uniform priors. The posterior probabilities of each of the the two outcomes given the other factor levels are calculated. See MacKay(2003, p. 460).
Bft2x2(X, div = 100, plotit = TRUE)
Bft2x2(X, div = 100, plotit = TRUE)
X |
a contingency table in the form of a 2 x 2 matrix with row and column names |
div |
optional: the number of divisions for the row and column variables for use in calculations (can be left at 100) |
plotit |
should plots be produced? (defaults to TRUE) |
Note that the rows of the 2 x 2 matrix are assumed to represent the "outcomes" and the columns the "treatments"—where these expressions are applicable. Note also that to obtain properly labeled plots the matrix has to be supplied with dimnames.
the matrix of div
x div
posterior probabilities that was plotted
Robert van Hulst
MacKay, D.J.C. 2003. Information Theory, Inference, and Learning Algorithms. Cambridge University Press, Cambridge.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(Glasses) Bft2x2(Glasses)
data(Glasses) Bft2x2(Glasses)
A total of eight models are fitted to a data set consisting of seven predictors. The response is the exact fit with a variable amount of zero-mean noise added. This is repeated a certain number of times (by default, 100 times). Plots of Bias^2 and variance vs. the number of parameters are produced.
BiasVarTO(times = 100)
BiasVarTO(times = 100)
times |
the number of repeats to average bias and variance over (default 100) |
none produced, the function produces two plots
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
This function just plots some Beta distributions with commonly used parameters
binPriorsPlot()
binPriorsPlot()
none produced, the function just produces one (compound) plot
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
binPriorsPlot()
binPriorsPlot()
These made-up data do respect the average clutch sizes (number of eggs laid in a single brood) and incubation periods that were observed in different European bird species with four different types of nests, as reported in Case(2000).
data(BirdsCS)
data(BirdsCS)
A data frame with 40 observations on the following 3 variables:
Nest
kind of nest, a factor with levels hole
,
roofed
, niche
, and open
Inc.Per
average duration of the incubation period (days)
ClutchSize
the typical number of eggs in a nest
Case, T.J. An Illustrated Guide to Theoretical Ecology. Oxford University Press, New York.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(BirdsCS) library(graphics) coplot(ClutchSize ~ Inc.Per | Nest, BirdsCS, panel=panel.smooth)
data(BirdsCS) library(graphics) coplot(ClutchSize ~ Inc.Per | Nest, BirdsCS, panel=panel.smooth)
Several exploratory plots are produced, after which this function calculates and plots the posterior densities of the treatment means and their differences. Pooled or separate variances can be specified. Note that this function uses Standard Improper Reference (SIR) priors.
BnNsir(formula, data, var.equal = TRUE, alpha = 0.05, plotit = TRUE, r = 10000)
BnNsir(formula, data, var.equal = TRUE, alpha = 0.05, plotit = TRUE, r = 10000)
formula |
the usual formula interface: response ~ factor |
data |
a data.frame containing the response and the factor variables |
var.equal |
should a pooled variance be used? Specify var.equal = FALSE if you want separate variances to be fitted |
alpha |
1 - level of credibility, so that for alpha = 0.05 (the default) credible intervals will have 95% credibility |
plotit |
are plots desired? |
r |
the number of samples of the posterior that should be taken |
none returned: the function is used for the plots and the printed information it produces
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(PlantGrowth) BnNsir(weight ~ group, PlantGrowth)
data(PlantGrowth) BnNsir(weight ~ group, PlantGrowth)
These data were collected by Mackowiak, Wasserman, and Levine(1992), and have been used, among others, by Ntzoufras(2009).
data(bodytemp)
data(bodytemp)
A data frame with 130 observations on the following 3 variables:
temperature
body temperature in degrees Fahrenheit
gender
a factor with levels 'female' and 'male'
heart.rate
heart rate in beats per minute
Mackowiak, P.A., Wasserman, S.S., and Levine, M.M.(1992) A critical appraisal of 98.6 degrees F, the upper limit of the normal body temperature, and other legacies of Carl Reinhold August Wunderlich. JASA 268, 1578–1580.
Ntzoufras, I.(2009) Bayesian Modeling Using Winbugs. Wiley, Hoboken, N.J.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(bodytemp) B2Nsir(temperature ~ gender, bodytemp)
data(bodytemp) B2Nsir(temperature ~ gender, bodytemp)
This function compares different linear models on the basis of their Bayes factors and by graphically comparing posterior model probabilities.
Bregbf(form.list, data, l=length(form.list))
Bregbf(form.list, data, l=length(form.list))
form.list |
a list of linear models, each expressed by a model formula, that should be compared; the models must all be applicable to the same data frame and use the same response variable |
data |
a data frame to be analyzed |
l |
the number of models to be compared; defaults to all models in the form.list |
Note that a list containing several appropriate models for the data frame should be prepared beforehand. See the example for how to do this.
A list with model parameter probabilities is silently returned.
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
## Not run: data(PlantGrowth) frmlst <- list( model0 = formula(weight ~ 1), model1 = formula(weight ~ group) ) Bregbf(form.list=frmlst, data=PlantGrowth) data(fev) frmlst.fev <- list( formula(FEV ~ Age), formula(FEV ~ Smoke), formula(FEV ~ Age + Smoke), formula(FEV ~ Age * Smoke) ) Bregbf(frmlst.fev, fev) ## End(Not run)
## Not run: data(PlantGrowth) frmlst <- list( model0 = formula(weight ~ 1), model1 = formula(weight ~ group) ) Bregbf(form.list=frmlst, data=PlantGrowth) data(fev) frmlst.fev <- list( formula(FEV ~ Age), formula(FEV ~ Smoke), formula(FEV ~ Age + Smoke), formula(FEV ~ Age * Smoke) ) Bregbf(frmlst.fev, fev) ## End(Not run)
The Bayesian “t-test” developed by Bernardo and Perez (2007) that calculates the Bayes-factor against the null hypothesis of no difference.
Bt.test(formula, data, plotit = TRUE)
Bt.test(formula, data, plotit = TRUE)
formula |
the usual formula interface: response ~ factor |
data |
a data.frame with the response values and the factor values for all samples; the factor can only have two factor levels |
plotit |
is plotted output required? |
none supplied: the function is used for the plotted and printed output it produces
Robert van Hulst
J. Bernardo and S. Perez. Comparing normal means: New methods for an old problem. Bayesian Analysis, 2:45–58, 2007.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(bodytemp) Bt.test(temperature ~ gender, bodytemp) Bt.test(heart.rate ~ gender, bodytemp)
data(bodytemp) Bt.test(temperature ~ gender, bodytemp) Bt.test(heart.rate ~ gender, bodytemp)
Batches of twenty larvae were exposed to increasing doses of insecticide, and the number of survivors and their sexes were noted. These data were reported by Collett(1991) and used by Venables and Ripley(1994 and later editions). They resulted from an experiment to study the toxicity of a pyrethroid insecticide to the tobacco budworm Heliothis virescens of different doses of the insecticide.
data(budworm)
data(budworm)
A data frame with 12 observations on the following 3 variables:
ldose
the log of the dose of the insecticide
dead
the number of budworms that were dead a day later
sex
a factor with two levels: “F” and “M”
Collett, D. 1991. Modelling Binary Data. Chapman and Hall, London.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
Venables, W.N. and Ripley, B.D. 1994. Modern Applied Statistics with S-PLUS. Springer Verlag, New York.
data(budworm) fit <- glm(cbind(dead, 20 - dead) ~ ldose, data=budworm, family=binomial) summary(fit)
data(budworm) fit <- glm(cbind(dead, 20 - dead) ~ ldose, data=budworm, family=binomial) summary(fit)
These made-up data illustrate the discrete form (contingency table form) of Simpson's paradox.
data(Clin)
data(Clin)
A three-dimensional array of frequencies with:
rows indicating "outcome" (either "death" or "cured"),
columns indicating "male" (either "Yes" or "No"), and
layers indicating "clinic" (either "A" or "B").
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(Clin) Clin[1,,] prop.table(Clin[1,,], 2)
data(Clin) Clin[1,,] prop.table(Clin[1,,], 2)
Nespolo et al.(2003) collected data on the metabolic rates (as measured by oxygen consumption) of crickets kept and acclimated at three different temperatures. Since the original data were not available and only a statistical summary was published, we simulated these data to approximately agree with the statistical summary.
data(crickets)
data(crickets)
A data frame with 292 observations on the following 3 variables:
VO2
oxygen consumption in l/h (a measure of basal metabolic rate)
mass
weight of the cricket in mg
temp
temperature in degrees C.
Nespolo et al., 2003.
Nespolo, R.F., Lardies, M.A., and Bozinovic, F. 2003. Intrapopulational variation in the standard metabolic rate of insects: repeatability, thermal dependence and sensitivity of (Q[10]) on oxygen consumption in a cricket. Journal of Experimental Biology 206, 4309–4315.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(crickets) crickets7 <- subset(crickets, crickets$temp==7) with(crickets7, scatter.smooth(mass, VO2))
data(crickets) crickets7 <- subset(crickets, crickets$temp==7) with(crickets7, scatter.smooth(mass, VO2))
An n x n contingency table is analyzed in frequentist, information-theoretical, likelihood, and Bayesian ways. Note that for the Bayesian analysis package LearnBayes needs to be installed.
CTA(X, extBayes = FALSE)
CTA(X, extBayes = FALSE)
X |
a matrix with non-negative integers representing the counts for the row-column levels |
extBayes |
should a Bayesian analysis with a near-independence prior (instead of only an independence prior) be done as well? Defaults to FALSE. |
none provided: the function is run for its graphical and numerical output
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(Smoking) CTA(Smoking)
data(Smoking) CTA(Smoking)
Charles Darwin(1876) provided data on the difference in the heights attained by selfed and crossed mother plants.
data(darwin)
data(darwin)
A data frame with 15 observations on the following variable:
difference
the difference in height in inches between each paired pair of offspring of a selfed and a crossed mother plant
Darwin, C.R. 1876. The effects of cross and self fertilisation in the vegetable kingdom. John Murray, London.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(darwin) with(darwin, qqnorm(difference) ) with(darwin, qqline(difference) )
data(darwin) with(darwin, qqnorm(difference) ) with(darwin, qqline(difference) )
Data from Johnson (1996) on human body fat: determined by under-water weight and several covariates to estimate it statistically.
data(Fat)
data(Fat)
A data frame with 252 observations on the following 19 variables:
Case
case number
PBF.B
percentage body fat estimated using Brozek's equation
PBF.S
percentage body fat estimated using Siri's equation
Dens
Density (gm/cm^3)
Age
Age (yrs)
Weight
Weight (lbs)
Height
Height (inches)
AI
Adiposity index = Weight/Height^2 (kg/m^2)
FFWt
Fat Free Weight using Brozek's formula (lbs)
Neck
Neck circumference (cm)
Chest
Chest circumference (cm)
Abd
Abdomen circumference (cm)
Hip
Hip circumference (cm)
Thigh
Thigh circumference (cm)
Knee
Knee circumference (cm)
Ankle
Ankle circumference (cm)
Biceps
Extended biceps circumference (cm)
FArm
Forearm circumference (cm)
Wrist
Wrist circumference (cm)
Johnson, R. 1996. Fitting percentage of body fat to simple body measurements. Journal of Statistics Education 2(1), 1–6.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(Fat) qqnorm(Fat$Height) qqline(Fat$Height)
data(Fat) qqnorm(Fat$Height) qqline(Fat$Height)
These data come from Rosner (2006), and represent forced expiratory volume (FEV) in l/s and several covariates.
data(fev)
data(fev)
A data frame with 654 observations on the following 6 variables:
Id
an identification code
Age
age in years
FEV
forced expiratory volume in l/s
Hgt
height in inches
Sex
gender: 0 for female, 1 for male
Smoke
smokes (1) or not (0)
Rosner, B. 2006. Fundamentals of Biostatistics. 6th ed. Duxbury Press.
data(fev) splom(fev[c(3, 2, 4, 5, 6)], main="fev data")
data(fev) splom(fev[c(3, 2, 4, 5, 6)], main="fev data")
Data from Heidelberger and Holland(2004) categorizing a random sample of 16 British juveniles on the basis of whether they were juvenile delinquents or not, and whether wore glasses or not.
data(Glasses)
data(Glasses)
A matrix with 16 counts cross-classified on Juvenile delinquency (rows) and the wearing of glasses (columns).
Heiberger, R.M. and Holland, B.(2004) Statistical Analysis and Data Display: An Intermediate Course with Examples in S-PLUS, R, and SAS. Springer, New York.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(Glasses) Bft2x2(Glasses)
data(Glasses) Bft2x2(Glasses)
function used to produce a Bayesian credible interval of a unimodal distribution of empirical values using the Highest Posterior Probability approach
HPDcrd(x, alpha = 0.05)
HPDcrd(x, alpha = 0.05)
x |
a vector of empirical values |
alpha |
1 - alpha is the desired level of credibility |
a vector of the lower and upper limits of the 95% credible interval calculated using a standard algorithm
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
HPDcrd(rnorm(1000))
HPDcrd(rnorm(1000))
Data on horseshoe crab morphology collected by Brockman(1996) and used by Agresti(2012).
data(HSCrab)
data(HSCrab)
A data frame with 173 observations on the following 5 variables:
Col
an indicator variable for the carapace color
spineW
coded width of the spine
Width
maximal width of the carapace (cm)
Satell
number of satellite males
Weight
weight in g
Brockman, H.J.(1996) Satellite male groups in horseshoe crabs, Limulus polyphemus Ethology 102(1), 1–21.
Agresti, A.(2012) Categorical Data Analysis (3rd ed.) Wiley, New York.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(HSCrab) plot(Weight ~ Width, col = Col, data = HSCrab)
data(HSCrab) plot(Weight ~ Width, col = Col, data = HSCrab)
When given the number of “successes” and the sample size this
function plots the normed likelihood of values of the binomial
parameter and calculates the likelihood ratio for a hypothesized
value and the maximum likelihood value for the sample, as well as an
approximate frequentist p-value.
L1prop(x, n, p.hypoth, pLset=0.05)
L1prop(x, n, p.hypoth, pLset=0.05)
x |
the number of sampling units with the feature |
n |
the number of sampling units examined |
p.hypoth |
the hypothesized probability |
pLset |
the desired likelihood for the likelihood interval |
none, the normed likelihood for different values of the binomial probability is plotted with the likelihood interval, and some information is printed
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
Pawitan, Y. 2001. In All Likelihood. Oxford University Press, Oxford.
L1prop(13, 78, 0.02)
L1prop(13, 78, 0.02)
When given the numbers of “successes” and the sample sizes for the two samples, this function plots the normed likelihoods of the two samples and calculates the likelihood ratio for two different models, one fitting two binomial parameters, and one fitting only one.
L2prop(x, n)
L2prop(x, n)
x |
a vector containing the 2 numbers of sampling units with the feature |
n |
a vector containing the 2 numbers of sampling units examined |
none, the inferred difference between the probabilities and its 95% credible interval are calculated and a plot is produced
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
L2prop(c(13, 22), c(78, 92))
L2prop(c(13, 22), c(78, 92))
Simon Newcom's measured in the late 1900's the time it took light to cover a certain distance. The data are reported in Stigler(1977) and have been widely used since to illustrate statistical inference.
data(lightspeed)
data(lightspeed)
A vector with 66 observations of the travel time of light.
Stigler, S.M. (1977) Do robust estimators work with real data? Annals of Statistics 5, 1055–1098.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(lightspeed) qqnorm(lightspeed) qqline(lightspeed)
data(lightspeed) qqnorm(lightspeed) qqline(lightspeed)
The LOOIC-value (like the non-Bayesian AIC-value) is a useful measure of model performance for model prediction.
looicplot(looiclist, modnames, perc = 90)
looicplot(looiclist, modnames, perc = 90)
looiclist |
a list of character-valued names of rstanarm model objects |
modnames |
a character-valued vector of model names for each of the models |
perc |
the percentage credibility for the credible intervals (defaults to 90%) |
None provided, but a printed list of looic-values, their standard errors, and credible intervals, and a dot plot with the same information are produced.
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
## Not run: data(budworm) Mbudworm1 <- stan_glm(formula = cbind(dead, 20 - dead) ~ ldose, family = binomial, data = budworm, prior = student_t(df = 7), prior_intercept = student_t(df = 7)) Mbudworm2 <- stan_glm(formula = cbind(dead, 20 - dead) ~ ldose * sex, family = binomial, data = budworm, prior = student_t(df = 7), prior_intercept = student_t(df = 7)) Mbudworm3 <- stan_glm(formula = cbind(dead, 20 - dead) ~ ldose + sex, family = binomial, data = budworm, prior = student_t(df = 7), prior_intercept = student_t(df = 7)) looicplot(looiclist = list("Mbudworm1", "Mbudworm2", "Mbudworm3"), modnames = c("~ ldose", "~ ldose + sex", "~ ldose * sex") ) ## End(Not run)
## Not run: data(budworm) Mbudworm1 <- stan_glm(formula = cbind(dead, 20 - dead) ~ ldose, family = binomial, data = budworm, prior = student_t(df = 7), prior_intercept = student_t(df = 7)) Mbudworm2 <- stan_glm(formula = cbind(dead, 20 - dead) ~ ldose * sex, family = binomial, data = budworm, prior = student_t(df = 7), prior_intercept = student_t(df = 7)) Mbudworm3 <- stan_glm(formula = cbind(dead, 20 - dead) ~ ldose + sex, family = binomial, data = budworm, prior = student_t(df = 7), prior_intercept = student_t(df = 7)) looicplot(looiclist = list("Mbudworm1", "Mbudworm2", "Mbudworm3"), modnames = c("~ ldose", "~ ldose + sex", "~ ldose * sex") ) ## End(Not run)
A strip chart of the first argument grouped by the second argument is produced. This function is useful for looking at experimental data with a numeric response and a factorial predictor.
meansplot(y, grp)
meansplot(y, grp)
y |
a vector of observed values |
grp |
a factor of the same length as the observation vector indicating the treatment under which each observation was obtained |
none returned: the function is used for the plot it produces
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(PlantGrowth) with(PlantGrowth, meansplot(weight, group))
data(PlantGrowth) with(PlantGrowth, meansplot(weight, group))
If experimental data on the sensitivity and the specificity of a diagnostic test are available, and the prevalence of the the condition is known with its raw data, then this function estimates the posterior probability of having the condition, with its 95% credible interval.
MedDiagn(x0, n0, x1, n1, x2, n2, N = 10000, alpha = 0.05, pdf = FALSE)
MedDiagn(x0, n0, x1, n1, x2, n2, N = 10000, alpha = 0.05, pdf = FALSE)
x0 |
prevalence raw data: number of people with a certain condition |
n0 |
number of people examined for that condition |
x1 |
sensitivity data: number of people with the disease for whom this test was positive |
n1 |
total number of people in the sensitivity sample |
x2 |
specificity raw data: number of people who did not have the disease who tested negative |
n2 |
total number of people in the specificity sample |
N |
number of cases to be simulated (best left at 10000 or greater |
alpha |
credibility required (default 95%) |
pdf |
set this to TRUE only if you want to keep a pdf-file of the posterior probability plot |
none returned: a plot and printed information are produced
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
MedDiagn(105, 35000, 72, 80, 640, 800)
MedDiagn(105, 35000, 72, 80, 640, 800)
Normality plots can be hard to judge if one is not experienced. This function plots a Normality plot for the data surrounded by eight other Normality plots for samples with the same mean and standard deviation that were randomly generated. The eight plots provide an idea of the variability to be expected in Normally distributed data.
nineplot(x)
nineplot(x)
x |
a vector of observations to be examined for Normality |
none produced: the function is used for the plot it produces
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
nineplot(rt(100, 2))
nineplot(rt(100, 2))
The negative predictive value (NPV) of a diagnostic test is the probability that someone with a negative diagnostic test for a condition does not have the condition. The NPV can easily be calculated from the prevalence, the sensitivity, and the specificity, but this function automates the procedure.
NPV(sens, spec, prev)
NPV(sens, spec, prev)
sens |
the sensitivity of the test |
spec |
the specificity of the test |
prev |
the prevalence of the disease |
the negative predictive value
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
NPV(0.9, 0.8, 0.003)
NPV(0.9, 0.8, 0.003)
A large sample of Normal-distributed data with more than 10% of the observations further than 1.5 times the IQR from the median shows signs of overdispersion, as recommended in Gelman et al., 2014.
overdispersionCheck(x)
overdispersionCheck(x)
x |
an input vector of reals without missing values |
The function prints the approximate percentage of observations that are further from the median than would be expected in a normal distribution.
Robert van Hulst
Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., and Rubin, D.B. 2014. Bayesian Data Analysis. Third Ed.. CRC Press
overdispersionCheck(rt(100, 1))
overdispersionCheck(rt(100, 1))
This function computes the approximate lower bound to the Bayes factor of the null hypothesis against the alternative, assuming equal odds of the null and the alternatlve.
p2BF(p)
p2BF(p)
p |
the frequentist p-value (which has to be less than 1/e or 0.37) |
the approximate lower bound of the Bayes factor of the null hypothesis against the alternative
the p-value should be less than 1/e (= 0.37).
Robert van Hulst
Sellke, T., Bayarri, M.J., and Berger, J.O. 2001. Calibration of p Values for Testing Precise Hypotheses. Am. Statistician 55(1) pp 62–71.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
p2BF(p = 0.05)
p2BF(p = 0.05)
This function computes the approximate lower bound to the posterior probability of the null hypothesis assuming equal odds of the null and the alternative. See Sellke et al.(2001) for the derivation, and note that the posterior probability of the null hypothesis is what many incorrectly assume the p-value is measuring.
p2minpp(p)
p2minpp(p)
p |
the frequentist p-value (which has to be less than 1/e or 0.37) |
the approximate lower bound of the posterior probability of the null hypothesis
the p-value should be less than 1/e (0.37).
Robert van Hulst
Sellke, T., Bayarri, M.J., and Berger, J.O. 2001. Calibration of p Values for Testing Precise Hypotheses. Am. Statistician 55(1) pp 62–71.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
p2minpp(p=0.05)
p2minpp(p=0.05)
The positive predictive value (PPV) of a dianostic test is the probability that someone with a positive diagnostic test for a condition does have the condition. The PPV can easily be calculated from the prevalence, the sensitivity, and the specificity, but this function automates the procedure.
PPV(sens, spec, prev)
PPV(sens, spec, prev)
sens |
the sensitivity of the test |
spec |
the specificity of the test |
prev |
the prevalence of the disease |
the positive predictive value of the test
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
PPV(0.9, 0.8, 0.003)
PPV(0.9, 0.8, 0.003)
Up to forty eight rats were exposed to the carcinogen retinyl acetate or to a placebo in their diet, after which the number of tumors they developed was evaluated.
data(rats)
data(rats)
A data frame with 71 observations on the following 2 variables:
y
number of rats that developed tumors
N
number of rats in group
Gail, M.H., Santner, T.J., and Brown, C.C. (1980) An analysis of comparative carcinogenesis experiments based on multiple times to tumor. Biometrics 36, 255–266.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
Given a critical value alpha, this function performs a Fisherian significance test of the null hypothesis at level p, reports the result of the test, as well as the lower and upper values of the corresponding confidence interval. See Kadane(2016) for the idea for this.
sigtestCI(p)
sigtestCI(p)
p |
the desired significance level |
Note that this function does not require any data: if a rare (as long as p is sufficiently small) event occurs, H[0] is deemed to be implausible, and rejected. If such an event does not occur, we can simply try to do the experiment again. A Neyman-Pearson hypothesis test does require data and also an alternative hypothesis. For a NP hypothesis test we can (and should) consider the power of the test (the probability of rejecting H[0] when H[a] is true).
A message informing the user if H0 was rejected or not and the lower and upper boundaries of the corresponding confidence interval.
Robert van Hulst
Kadane, J.B. 2016. Beyond hypothesis testing. Entropy 18, 199.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
sigtestCI(p=0.05)
sigtestCI(p=0.05)
The data are from a retrospective study that compared mortality due to a heart infarct in people who smoked and sex-matched controls who did not.
data(Smoking)
data(Smoking)
A matrix with 781 observations cross-classified on the following 2 factors: “Infarct” (”Yes” or “Control”, rows), and “EverSmoked” (”Yes” or “No”, columns).
unknown
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
This function converts the successes and totals vectors required as input for function B2props to a 2x2 contingency table for input to CTA or Bft2x2.
sn2ft2x2(s, n)
sn2ft2x2(s, n)
s |
a vector of length 2 of successes |
n |
a vector of length 2 of numbers of trials |
a 2 x 2 contingency table equivalent to the two arguments
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
sn2ft2x2(c(47, 59), c(120, 125))
sn2ft2x2(c(47, 59), c(120, 125))
A total of 433 persons were tested for hypertension and checked for whether they were smokers, obese, or snored. The data are in Altman(1991).
data(Snoring)
data(Snoring)
A data frame with 8 observations on the following 5 variables:
smoking
did the person smoke (1) or not (0)?
obese
was the person obese (1) or not (0)?
snoring
did the person snore (1) or not (0)?
n
the number of persons observed with these covariates
hypert
did the person suffer from hypertension (1) or not (0)?
Altman, D.G. 1991. Practical Statistics for Medical Research. Chapman \& Hall, London.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(Snoring) fit <- glm(cbind(hypert, n - hypert) ~ smoking + obese + snoring, family=binomial, data=Snoring) summary(fit)
data(Snoring) fit <- glm(cbind(hypert, n - hypert) ~ smoking + obese + snoring, family=binomial, data=Snoring) summary(fit)
These data came from a designed experiment reported in Sokal and Rohlf(1995), box 9.4. The growth (in arbitrary units) of pea sections grown in tissue culture on five different sugars was replicated ten times.
data(SRb94)
data(SRb94)
A data frame with 50 observations on the following 2 variables:
L
length difference in mm
Treatm
a factor with levels "Contr", "fruct.", "gluc.", "gluc&fruct.", and "sucr."
Sokal, R.R., and Rohlf, F.J. Biometry. Freeman, New York.
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
data(SRb94) with(SRb94, meansplot(L, Treatm))
data(SRb94) with(SRb94, meansplot(L, Treatm))
The sum of squares of the input vector is returned.
SSQ(x)
SSQ(x)
x |
a vector of numbers without missing values |
the sum of squares of x
Robert van Hulst
SSQ(x = rnorm(n=100))
SSQ(x = rnorm(n=100))
Produces a dotchart with error bars as summary of a dataframe with model names (‘modnames’), LOOIC-values (‘looic’), standard errors (‘se’), lower values (‘lwr’), and upper values (‘upr’) .
sumchart(df, rownames, groups, perc)
sumchart(df, rownames, groups, perc)
df |
data.frame name |
rownames |
model names |
groups |
row names |
perc |
the percentage of credibility desired |
A plot is produced.
Robert van Hulst
van Hulst, R. 2018. Evaluating Scientific Evidence. ms.
Rats were fed diets with different quantities of protein from either animal or plant sources. The weight gained at the end of the experiment was the response variable.
data("weightgain")
data("weightgain")
A data frame with 40 observations on the following 3 variables
source
source of protein given, a factor with levels
Beef
and Cereal
type
amount of protein given, a factor with levels High
and Low
weightgain
weight gain in grams
Hand, D.J., Daly, F., Lunn, A.D., McConway, K.J. and Ostrowski, E. 1994. A Handbook of Small Datasets, Chapman and Hall, London.
data("weightgain") with(weightgain, table(source, type))
data("weightgain") with(weightgain, table(source, type))