| Title: | Likelihood Analyses for Common Statistical Tests |
|---|---|
| Description: | A collection of functions that calculate the log likelihood (support) for a range of statistical tests. Where possible the likelihood function and likelihood interval for the observed data are displayed. The evidential approach used here is based on the book "Likelihood" by A.W.F. Edwards (1992, ISBN-13 : 978-0801844430), "Statistical Evidence" by R. Royall (1997, ISBN-13 : 978-0412044113), S.N. Goodman & R. Royall (2011) <doi:10.2105/AJPH.78.12.1568>, "Understanding Psychology as a Science" by Z. Dienes (2008, ISBN-13 : 978-0230542310), S. Glover & P. Dixon <doi:10.3758/BF03196706> and others. This package accompanies "Evidence-Based Statistics" by P. Cahusac (2020, ISBN-13 : 978-1119549802) <doi:10.1002/9781119549833>. |
| Authors: | Peter Cahusac [aut, cre]
|
| Maintainer: | Peter Cahusac <[email protected]> |
| License: | GPL-2 |
| Version: | 1.1.5 |
| Built: | 2024-12-22 06:20:51 UTC |
| Source: | CRAN |
This function calculates supports for independent samples ANOVA. One support is for the model of group means against the null (no grouping), for the first contrast versus the group means model, and the other for 2 contrasts. Both contrasts should be either NULL or specified. If the contrasts use the default of NULL, then it calculates a linear versus a quadratic contrast. The corrected support is given for groups versus null, using Akaike's correction (Hurvich & Tsai (1989)). No correction is necessary for the two contrasts' support since they both involve 1 parameter. Conventional frequentist F and p value statistics are given for the overall analysis and for contrast 1. Unequal group sizes are accommodated.
L_1way_ANOVA(data, group, contrast1=NULL, contrast2=NULL, verb=TRUE)L_1way_ANOVA(data, group, contrast1=NULL, contrast2=NULL, verb=TRUE)
data |
a (non-empty) numeric vector of data values. |
group |
an integer vector the same length as data, coding for k groups. |
contrast1 |
first contrast, default = NULL. |
contrast2 |
second contrast, default = NULL. |
verb |
show output, default = TRUE. |
$S.12c - corrected support for groups hypothesis versus null.
$S.12 - uncorrected support for groups hypothesis versus null.
S.1mc - support for contrast 1 versus the group means model.
$S.1m = uncorrected support for contrast 1 vs group means.
$S.cont.12 - support for contrast 1 versus contrast 2.
$contrast1 - first contrast.
$contrast2 - second contrast.
$gp.means - group means.
$df - degrees of freedom for groups and error.
$F.val - F value for overall ANOVA analysis.
$P.val - p value for overall analysis.
$eta.sq - eta-squared.
$Fval.c1 - F value for contrast 1.
$df.1 - degrees of freedom for contrast.
$P.val1 - p value for contrast 1.
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Hurvich CM, Tsai C-L. Regression and time series model selection in small samples. Biometrika. 1989; 76(2):297.
Dixon P. The effective number of parameters in post hoc models. Behavior Research Methods. 2013; 45(3):604.
Dixon P. The p-value fallacy and how to avoid it. Canadian Journal of Experimental Psychology/Revue canadienne de psychologie expérimentale. 2003; 57(3):189.
Glover S, Dixon P. Likelihood ratios: a simple and flexible statistic for empirical psychologists. Psychonomic Bulletin and Review. 2004; 11(5):791.
# fitness example, p 81 dat <- c(7, 5, 9, 8, 3, 12, 10, 8, 7, 9, 5, 7, 7, 6, 4, 8, 12, 9, 7, 8, 3, 2, 7, 6, 8, 6, 5, 3, 4, 3, 4, 3, 3, 1, 2, 5, 7, 6, 8, 7) gp <- as.factor(rep(1:4,each=10)) gp = gl(4,10,40, labels=c("6 hr", "3 hr", "Sports club", "Video games")) contrast1 <- c(-3, -1, 1, 3) # linear contrast2 <- c(1, -1, -1, 1) # quadratic contrast3 <- c(1, 1, -1, -1) L_1way_ANOVA(dat,gp,contrast3, contrast1)# fitness example, p 81 dat <- c(7, 5, 9, 8, 3, 12, 10, 8, 7, 9, 5, 7, 7, 6, 4, 8, 12, 9, 7, 8, 3, 2, 7, 6, 8, 6, 5, 3, 4, 3, 4, 3, 3, 1, 2, 5, 7, 6, 8, 7) gp <- as.factor(rep(1:4,each=10)) gp = gl(4,10,40, labels=c("6 hr", "3 hr", "Sports club", "Video games")) contrast1 <- c(-3, -1, 1, 3) # linear contrast2 <- c(1, -1, -1, 1) # quadratic contrast3 <- c(1, 1, -1, -1) L_1way_ANOVA(dat,gp,contrast3, contrast1)
This function calculates the support for one-way categorical data (multinomial), also gives chi-squared and likelihood ratio test (G) statistics. If there are only 2 categories then binomial information is given too with likelihood interval, including the likelihood-based % confidence interval. Support for the variance being more different than expected (Edwards p 187, Cahusac p 158) is also calculated. It uses the optimize function to locate desired limits for both intervals.
L_1way_cat(obs, exp.p=NULL, L.int=2, alpha=0.05, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)L_1way_cat(obs, exp.p=NULL, L.int=2, alpha=0.05, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)
obs |
a vector containing the number of counts in each category. |
exp.p |
a vector containing expected probabilities. If NULL then this is 1/#cats. |
L.int |
likelihood interval given as support values, e.g. 2 or 3, default = 2. |
alpha |
the significance level used, 1 - alpha interval calculated, default = 0.05. |
toler |
the desired accuracy using optimise, default = 0.0001. |
logplot |
plot vertical axis as log likelihood, default = FALSE |
supplot |
set minimum likelihood display value in plot, default = -10 |
verb |
show output, default = TRUE. |
$S.val - support for one-way observed versus expected.
$uncorrected.sup - uncorrected support.
$df - degrees of freedom for table.
$observed - observed counts.
$exp.p - expected probabilities.
$too.good - support for the variance of counts being more different than expected.
$chi.sq - chi-squared value.
$p.value - p value for chi-squared.
$LR.test = the likelihood ratio test statistic.
$lrt.p = the p value for the likelihood ratio test statistic
Additional outputs for binomial:
$prob.val - MLE probability from data.
$succ.fail - number of successes and failures.
$like.int - likelihood interval.
$like.int.spec - specified likelihood interval in units of support.
$conf.int - likelihood-based confidence interval.
$alpha.spec - specified alpha for confidence interval.
$err.acc - error accuracy for optimize function.
Aitkin, M. et al (1989) Statistical Modelling in GLIM, Clarendon Press, ISBN : 978-0198522041
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Royall, R. M. (1997). Statistical evidence: A likelihood paradigm. London: Chapman & Hall, ISBN : 978-0412044113
Edwards, A.W.F. (1992) Likelihood, Johns Hopkins Press, ISBN : 978-0801844430
# example for binomial, p 123 obs <- c(6,4) L_1way_cat(obs, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb = TRUE) # example for multinomial, p 134 obs <- c(60,40,100) exp <- c(0.25,0.25,0.5) L_1way_cat(obs, exp.p=exp, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb = TRUE)# example for binomial, p 123 obs <- c(6,4) L_1way_cat(obs, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb = TRUE) # example for multinomial, p 134 obs <- c(60,40,100) exp <- c(0.25,0.25,0.5) L_1way_cat(obs, exp.p=exp, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb = TRUE)
This function calculates support for the treatment means versus the null model, Type of correction for the parameters in the null versus the means model can be selected. support for the first contrast versus the group means model and the support for the first versus the second contrast. Also gives the F, p and partial eta-squared values for the overall analysis, and F and p for the first contrast. Assumes sphericity.
L_1way_RM_ANOVA(dat, group, ID, correct=1, contrast1=NULL, contrast2=NULL, verb=TRUE)L_1way_RM_ANOVA(dat, group, ID, correct=1, contrast1=NULL, contrast2=NULL, verb=TRUE)
dat |
a (non-empty) numeric vector of data values. |
group |
a vector the same length as data, coding for k groups. |
ID |
is an identifier for each case. |
correct |
specifies the correction: 0 = none, 1 = Occam's bonus (default), 2 = AIC |
contrast1 |
first contrast, default = NULL. |
contrast2 |
second contrast, default = NULL. |
verb |
show output, default = TRUE. |
$S.12 - support for groups means hypothesis versus null, uncorrected
$S.12c - corrected
$capplied - type of correction
$S.1m = support for first contrast versus means model.
$S.cont.12 - support for first versus second contrast.
$contrast1 - first contrast.
$contrast2 - second contrast.
$gp.means - group means.
$df - degrees of freedom for ANOVA.
$F.val - F value for overall analysis.
$P.val - p value for the overall analysis.
Fval.c1 - F value for the first contrast.
$df.1 - degrees of freedom for the contrast.
$P.val1 - p value for the contrast.
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Dixon P. The effective number of parameters in post hoc models. Behavior Research Methods. 2013; 45(3):604.
Dixon P. The p-value fallacy and how to avoid it. Canadian Journal of Experimental Psychology/Revue canadienne de psychologie expérimentale. 2003; 57(3):189.
Glover S, Dixon P. Likelihood ratios: a simple and flexible statistic for empirical psychologists. Psychonomic Bulletin and Review. 2004; 11(5):791.
# sleep data example, p 97 dat <- c(0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0, 2, 1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4, 1.5, 1.4, 0.0, -0.7, 0.5, 5.1, 5.7, 1.5, 4.7, 3.5) treat <- gl(3,10,30) patients <- gl(10,1,30) contrast1 <- c(-1, 0, 1) # linear contrast2 <- c(-2, 1, 1) # 1st vs treatments 2 & 3 m=L_1way_RM_ANOVA(dat, treat, patients, correct=1, contrast1, contrast2) m# sleep data example, p 97 dat <- c(0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0, 2, 1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4, 1.5, 1.4, 0.0, -0.7, 0.5, 5.1, 5.7, 1.5, 4.7, 3.5) treat <- gl(3,10,30) patients <- gl(10,1,30) contrast1 <- c(-1, 0, 1) # linear contrast2 <- c(-2, 1, 1) # 1st vs treatments 2 & 3 m=L_1way_RM_ANOVA(dat, treat, patients, correct=1, contrast1, contrast2) m
This function calculates several different supports for independent samples. Effect size (Cohen's d) and a second alternative hypothesis value can be specified. The maximum support is the support for the observed mean versus the null value. The support for the specified d versus the null is also calculated. If a second hypothesis value is specified (in units of the original measurements) then two further supports are calculated: d versus 2nd alternative hypothesis, and 2nd alternative hypothesis versus the null. The likelihood curve graphic with MLE and specified hypothesis values is produced. Finally, the requested likelihood interval is provided. The t, p and observed d values for the test against the null are given. If variances are specified as unequal then uses Welch's test where homogeneity of variance is not required.
L_2S_ttest(data, group, veq=0, null=0, d=0.5, alt.2=NULL, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)L_2S_ttest(data, group, veq=0, null=0, d=0.5, alt.2=NULL, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)
data |
a (non-empty) numeric vector of data values. |
group |
an integer vector the same length as data, coding for 2 groups. |
veq |
whether variances are equal: 1 = Yes, 0 = No, default = 0. |
null |
value for the null hypothesis, default = 0. |
d |
Cohen's effect size, default = 0.5. |
alt.2 |
value for an alternative hypothesis, in units used for data, default = NULL. |
L.int |
likelihood interval given as support values, e.g. 2 or 3, default = 2. |
toler |
the desired accuracy using optimise, default = 0.0001. |
logplot |
plot vertical axis as log likelihood, default = FALSE |
supplot |
set minimum likelihood display value in plot, default = -10 |
verb |
show output, default = TRUE. |
$obs.diff - the observed difference in means.
$df - degrees of freedom.
$var.eq - if not equal (0) then Welch's test used.
$alt.H1 - mean value according to specified d.
$alt.H2 - specified second hypothesis value.
$S_max - maximum support for observed mean difference against the null.
$S_10 - support for d versus null.
$S_12 - support for d versus specified second hypothesis.
$S_20 - support for second hypothesis versus the null.
$like.int - likelihood interval.
$L.int.spec - specified likelihood interval in units of support.
$null.value - null value
$t.val - t value for test against null.
$p.val - p value for test against null.
$d.obs - observed effect size (from null).
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Baguley, T. (2012) Serious Stats, Palgrave Macmillan, ISBN: 978-0230577183
Edwards, A.W.F. (1992) Likelihood, Johns Hopkins Press, ISBN : 978-0801844430
Royall, R. M. (1997) Statistical Evidence: A Likelihood Paradigm, Chapman & Hall, ISBN : 978-0412044113
Royall, R. M. (2000). On the probability of observing misleading statistical evidence. Journal of the American Statistical Association, 95, 760.
# using a variation on Gosset's original additional hours of sleep data, p 59 mysample <- c(0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0) treat <- rep(1:0,each=5) L_2S_ttest(mysample, treat, veq=0, null=0, d=0.5, alt.2=2, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)# using a variation on Gosset's original additional hours of sleep data, p 59 mysample <- c(0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0) treat <- rep(1:0,each=5) L_2S_ttest(mysample, treat, veq=0, null=0, d=0.5, alt.2=2, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)
This function calculates supports for two-way categorical data. This consists of the support for the interaction and the two main effects. Support for the interaction being closer or worse (different variance) than expected (Edwards p 187, Cahusac p 158) is calculated. The support for trend across the columns is given (assuming the levels for columns are ordered), and conventional p value for trend. Finally, Chi-squared and likelihood ratio test (G) statistics are given.
L_2way_cat(table, verb=TRUE)L_2way_cat(table, verb=TRUE)
table |
a 2 x 2 matrix or contingency table containing counts. |
verb |
show output, default = TRUE. |
$S.int - support for the interaction.
$df - the degrees of freedom for the interaction.
$S.int.unc - the uncorrected support for the interaction.
$S.Main.rows - support for the rows main effect.
$S.Main.cols - support for the columns main effect.
$S.Mr.uncorr - uncorrected support for rows main effect.
$S.Mc.uncorr - uncorrected support for the columns main effect.
$df.rows - degrees of freedom for rows.
$df.cols - degrees of freedom for columns.
$S.total - support for the whole table.
$S.trend - support for the trend across columns (if ordered).
$too.good - support for the variance being different from expected.
$observed - the observed table frequencies.
$expected - the expected values for null hypothesis of no interaction.
$residuals - the Pearson residuals.
$LR.test = the likelihood ratio test statistic.
$lrt.p - the p value for likelihood ratio test.
$chi.sq - chi-squared value.
$p.value - p value for chi-squared.
$trend.p - p value for trend (from chi-squared dist.).
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Royall, R. M. (1997). Statistical evidence: A likelihood paradigm. London: Chapman & Hall, ISBN : 978-0412044113
Edwards, A.W.F. (1992) Likelihood, Johns Hopkins Press, ISBN : 978-0801844430
# S. mansoni eggs in stools example, p 151 eggs <- as.table(rbind(c(14, 16, 14, 7, 6), c(87, 33, 66, 34, 11))) dimnames(eggs) = list("Infested" = c("Positive","Negative"), "Age Group" = c("0-","10-", "20-", "30-", "40-")) L_2way_cat(eggs) # or as a matrix eggs <- as.matrix(c(14, 87, 16, 33, 14, 66, 7, 34, 6, 11)) dim(eggs) <- c(2,5) L_2way_cat(eggs)# S. mansoni eggs in stools example, p 151 eggs <- as.table(rbind(c(14, 16, 14, 7, 6), c(87, 33, 66, 34, 11))) dimnames(eggs) = list("Infested" = c("Positive","Negative"), "Age Group" = c("0-","10-", "20-", "30-", "40-")) L_2way_cat(eggs) # or as a matrix eggs <- as.matrix(c(14, 87, 16, 33, 14, 66, 7, 34, 6, 11)) dim(eggs) <- c(2,5) L_2way_cat(eggs)
This function calculates supports for independent samples ANOVA. One support is for the full model versus null (no factors), and the second is for full model versus main effects. Two contrasts can be specified which can be used to explore interactions. Each should be given as a vector arranged as means for factor1 changing first (see example). If only the first contrast is specified then this is compared to the main effects model. If a second contrast is specified then the first contrast is compared to it. Corrected support is given where appropriate, using Akaike's correction (Hurvich & Tsai (1989)). No correction is necessary for the two contrasts support since they both involve 1 parameter. Unequal group sizes are accommodated, using type III sums of squares. F, p and partial eta-squared values are given for the two factors and their interaction.
L_2way_Factorial_ANOVA(data, factor1, factor2, contrast1=NULL, contrast2=NULL, verb=TRUE)L_2way_Factorial_ANOVA(data, factor1, factor2, contrast1=NULL, contrast2=NULL, verb=TRUE)
data |
a (non-empty) numeric vector of data values. |
factor1 |
a vector the same length as data, coding the first factor. |
factor2 |
a vector the same length as data, coding the second factor. |
contrast1 |
first contrast, with values for factor1 changing first, default = NULL. |
contrast2 |
second contrast, default = NULL. |
verb |
show output, default = TRUE. |
$S.12c - corrected support for full model versus null.
$S.12 - uncorrected support for full model versus null.
$S_FMc - corrected support for full model versus main effects model.
$S.FM - uncorrected support for full versus main effects.
$S.c1.Mc - corrected support for first contrast versus main effects model.
$S.c1.M - uncorrected support for first contrast versus main effects.
$S.c1.c2 - support for first versus second contrast.
$Means - 2 way table of means.
$df - degrees of freedom for the ANOVA.
$F.f1 - F value for first factor main effect.
$Pval.f1 - P value for first factor main effect.
$eta.sq.1 - partial eta-squared for first factor main effect.
$F.f2 - F value for second factor main effect.
$Pval.f2 - P value for second factor main effect.
$eta.sq.2 - partial eta-squared for second factor main effect.
$F.int - F value for interaction.
$Pval.int - P value for interaction.
$eta.sq.12 - partial eta-squared for the interaction.
$F.val.c1 - F value for first contrast.
$P.val.c1 - P value for first contrast.
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Hurvich CM, Tsai C-L. Regression and time series model selection in small samples. Biometrika. 1989; 76(2):297.
Dixon P. The effective number of parameters in post hoc models. Behavior Research Methods. 2013; 45(3):604.
Dixon P. The p-value fallacy and how to avoid it. Canadian Journal of Experimental Psychology/Revue canadienne de psychologie expérimentale. 2003; 57(3):189.
Glover S, Dixon P. Likelihood ratios: a simple and flexible statistic for empirical psychologists. Psychonomic Bulletin and Review. 2004; 11(5):791.
# blood clotting times example, p 91 time <- c(6.4, 4.6, 6.4, 5.6, 5.9, 6.1, 6.3, 4.5, 4.8, 6.6, 7, 9.3, 7.9, 9.4, 8.2, 4.4, 4.2, 5, 6.9, 4.5, 4, 4.3, 6.9, 5.5, 5.8, 4.4, 4.2, 5.1, 6.9, 4.5) Treatment = gl(3,5,30, labels=c("T1","T2","T3")) Health = gl(2,15,30, labels=c("Hemophiliac","Normal")) L_2way_Factorial_ANOVA(time, Treatment, Health) contrast1 <- c(-1, -1, 5, -1, -1, -1) # interaction Hemo T3 higher than others L_2way_Factorial_ANOVA(time, Treatment, Health, contrast1) contrast2 <- c(-1, -1, -1, 1, 1, 1) # main effect of health status (Hemo higher than Normal) m=L_2way_Factorial_ANOVA(time, Treatment, Health, contrast1, contrast2) m #show outputs# blood clotting times example, p 91 time <- c(6.4, 4.6, 6.4, 5.6, 5.9, 6.1, 6.3, 4.5, 4.8, 6.6, 7, 9.3, 7.9, 9.4, 8.2, 4.4, 4.2, 5, 6.9, 4.5, 4, 4.3, 6.9, 5.5, 5.8, 4.4, 4.2, 5.1, 6.9, 4.5) Treatment = gl(3,5,30, labels=c("T1","T2","T3")) Health = gl(2,15,30, labels=c("Hemophiliac","Normal")) L_2way_Factorial_ANOVA(time, Treatment, Health) contrast1 <- c(-1, -1, 5, -1, -1, -1) # interaction Hemo T3 higher than others L_2way_Factorial_ANOVA(time, Treatment, Health, contrast1) contrast2 <- c(-1, -1, -1, 1, 1, 1) # main effect of health status (Hemo higher than Normal) m=L_2way_Factorial_ANOVA(time, Treatment, Health, contrast1, contrast2) m #show outputs
This function calculates the support for a correlation from 2 vectors of data. An expected correlation can be specified and the support calculated for this relative to the observed and the null (which is assumed to be 0, but can also be specified) values. A likelihood function is plotted for the obtained correlation with a likelihood interval added and expected correlation, if specified. Conventional p value is also given.
L_corr(xv, yv, null=0, exp.r=NULL, L.int=2, alpha=.05, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)L_corr(xv, yv, null=0, exp.r=NULL, L.int=2, alpha=.05, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)
xv |
a numeric vector. |
yv |
a numeric vector the same length as xv. |
null |
the null value, default = 0. |
exp.r |
a specified correlation (could be expected value for the study), default = NULL. |
L.int |
likelihood interval given as support values, e.g. 2 or 3, default = 2. |
alpha |
the significance level used, 1 - alpha interval calculated, default = 0.05. |
toler |
the desired accuracy using optimise, default = 0.0001. |
logplot |
plot vertical axis as log likelihood, default = FALSE |
supplot |
set minimum likelihood display value in plot, default = -10 |
verb |
show output, default = TRUE. |
$obs.r - observed correlation.
$S.0 - support for observed correlation versus the null.
$S.1 - support for the specified correlation versus observed correlation.
$S.10 - support for the specified correlation versus the null value.
$exp.r - the specified correlation.
$N - the sample size.
$p.value - the p value for significance test versus 0.
$like.int - the likelihood interval.
$like.int.spec - the specified likelihood interval in terms of support.
$conf.int - the % confidence interval for the correlation.
$alpha.spec - the specified alpha for the % confidence interval.
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Royall, R. M. (1997). Statistical evidence: A likelihood paradigm. London: Chapman & Hall, ISBN : 978-0412044113
Edwards, A.W.F. (1992) Likelihood, Johns Hopkins Press, ISBN : 978-0801844430
# for heptathlon example, p 104 m200 <- c(22.6, 23.7, 23.1, 23.6, 23.6, 23.6, 25.5, 23.9, 24.5, 23.9, 24.9, 24.8, 24.7, 25.0, 24.6, 24.9, 25.0, 25.6, 24.8, 25.5, 25.7, 24.9, 26.6, 25.2, 26.2) m800 <- c(128.5, 126.1, 124.2, 132.5, 134.7, 132.5, 138.5, 127.9, 133.7, 132.2, 136.1, 142.8, 125.8, 131.5, 137.1, 134.9, 146.7, 133.9, 146.4, 144.0, 133.4, 138.0, 139.2, 137.3, 163.4) m=L_corr(m200, m800, null=0, exp.r=.4, L.int=3, alpha=.05, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE) m #Note: the support for observed vs 0 is different from book (5.776 vs 5.700) #due to differences in calculation of r by Excel and R# for heptathlon example, p 104 m200 <- c(22.6, 23.7, 23.1, 23.6, 23.6, 23.6, 25.5, 23.9, 24.5, 23.9, 24.9, 24.8, 24.7, 25.0, 24.6, 24.9, 25.0, 25.6, 24.8, 25.5, 25.7, 24.9, 26.6, 25.2, 26.2) m800 <- c(128.5, 126.1, 124.2, 132.5, 134.7, 132.5, 138.5, 127.9, 133.7, 132.2, 136.1, 142.8, 125.8, 131.5, 137.1, 134.9, 146.7, 133.9, 146.4, 144.0, 133.4, 138.0, 139.2, 137.3, 163.4) m=L_corr(m200, m800, null=0, exp.r=.4, L.int=3, alpha=.05, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE) m #Note: the support for observed vs 0 is different from book (5.776 vs 5.700) #due to differences in calculation of r by Excel and R
This function calculates the support for the efficacy, the likelihood interval and the likelihood-based confidence interval. It uses the optimize function to locate desired limits and their error.
L_efficacy(a, n, null=0, exp.eff=NULL, L.int=2, alpha=0.05, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)L_efficacy(a, n, null=0, exp.eff=NULL, L.int=2, alpha=0.05, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)
a |
the number of affected in control group. |
n |
total number of participants. |
null |
the null value for efficacy, if no effect then it would be 0, default = 0. |
exp.eff |
the expected or hypothesized efficacy, default = NULL. |
L.int |
likelihood interval given as support values, e.g. 2 or 3, default = 2. |
alpha |
the significance level used, 1 - alpha interval calculated, default = 0.05. |
toler |
the desired accuracy using optimise, default = 0.0001. |
logplot |
plot vertical axis as log likelihood, default = FALSE |
supplot |
set minimum likelihood display value in plot, default = -10 |
verb |
show output, default = TRUE. |
$S.val - support for the observed efficacy versus the null value.
$obs.eff - the observed efficacy.
$null - the null efficacy.
$exp.eff - expected efficacy as specified.
$S.exp.vsObs - support for expected efficacy versus observed.
$S.exp.versus.null - support for the expected efficacy versus the null.
$L.int - the likelihood interval for the observed efficacy.
$S_int - the specified likelihood interval.
$observed - observed numbers affected in control and intervention groups.
$expected - expected numbers according to the null.
$chi.sq - chi-squared statistic.
$p.value - p value associated with chi-squared statistic.
$df - degrees of freedom for chi-squared.
$residuals - the Pearson residuals.
$conf.int - likelihood-based confidence interval according to specified alpha.
$alpha - specified alpha for confidence interval.
$all.err.acc - error accuracy for each application of the optimize function.
Aitkin, M. et al (1989) Statistical Modelling in GLIM, Clarendon Press, ISBN : 978-0198522041
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Royall, R. M. (1997). Statistical evidence: A likelihood paradigm. London: Chapman & Hall, ISBN : 978-0412044113
Edwards, A.W.F. (1992) Likelihood, Johns Hopkins Press, ISBN : 978-0801844430
# pfizer covid-19 efficacy 2020 m = L_efficacy(a = 86, n = 94, null=0.8, exp.eff=0.95, L.int=2, alpha=0.05, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE) m# pfizer covid-19 efficacy 2020 m = L_efficacy(a = 86, n = 94, null=0.8, exp.eff=0.95, L.int=2, alpha=0.05, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE) m
This function calculates the supports for multiple logistic regression. A binary dependent variable is entered into the function, followed by up to 6 predictor variables (which need to be dummy coded if nominal and more than 2 levels). Outputs give the overall support for the full model versus the null (constant) model, supports for each of the predictor variables. Outputs include the usual chi-squared and p values.
L_logistic_regress(yv, p1, p2=NULL, p3=NULL, p4=NULL, p5=NULL, p6=NULL, verb=TRUE)L_logistic_regress(yv, p1, p2=NULL, p3=NULL, p4=NULL, p5=NULL, p6=NULL, verb=TRUE)
yv |
a binomial numeric vector for dependent variable. |
p1 |
vector for predictor variable, same length as yv. |
p2 |
vector for predictor variable, same length as yv, default = NULL. |
p3 |
vector for predictor variable, same length as yv, default = NULL. |
p4 |
vector for predictor variable, same length as yv, default = NULL. |
p5 |
vector for predictor variable, same length as yv, default = NULL. |
p6 |
vector for predictor variable, same length as yv, default = NULL. |
verb |
show output, default = TRUE. |
$S.FNc - support for the full versus null (constant value) model.
$S.each - support for each of the predictors, from first to last.
$chi.sq.FN - chi-squared for full versus null model.
$df - degrees of freedom for chi-squared.
$chi.sq.FN.p - p value for chi-squared.
$p.vals - p values for each of the predictors
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
1 Akaike H. A new look at the statistical model identification. IEEE Transactions on Automatic Control. 1974;19(6):716.
Glover S, Dixon P. Likelihood ratios: a simple and flexible statistic for empirical psychologists. Psychonomic Bulletin & Review. 2004;11(5):791.
Tabachnick BG, Fidell LS. Using Multivariate Statistics. Boston: Pearson Education; 2007.
# prescription errors example, p 114 p_error <- c(rep(1,6),rep(0,9)) score <- c(4, 5, 6, 5, 4, 6, 6, 4, 5, 8, 9, 7, 10, 8, 9) med1 <- c(1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1) med2 <- c(0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0) m1 = L_logistic_regress(p_error, score, med1, med2) m1# prescription errors example, p 114 p_error <- c(rep(1,6),rep(0,9)) score <- c(4, 5, 6, 5, 4, 6, 6, 4, 5, 8, 9, 7, 10, 8, 9) med1 <- c(1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1) med2 <- c(0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0) m1 = L_logistic_regress(p_error, score, med1, med2) m1
This function calculates the support for an OR from a 2 x 2 categorical data table. An expected OR can be specified and the support calculated for this relative to the observed and null (which is assumed to be 1, but can also be specified) values. A likelihood function is plotted for the obtained OR with a specified likelihood interval, and expected OR, if specified. The log likelihood plot can optionally be given instead. Chi-squared and likelihood ratio test (G) statistics are also provided and a likelihood-based % confidence interval. It uses the optimize function to locate desired limits for both intervals and other support calculations.
L_OR(table, null=1, exp.OR=NULL, L.int=2, alpha=0.05, cc=FALSE, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)L_OR(table, null=1, exp.OR=NULL, L.int=2, alpha=0.05, cc=FALSE, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)
table |
a 2 x 2 matrix or contingency table containing counts. |
null |
the value against which the obtained OR is tested, default = 1. |
exp.OR |
an expected or hypothetical OR. |
L.int |
likelihood interval given as support values, e.g. 2 or 3, default = 2. |
alpha |
the significance level used, 1 - alpha interval calculated, default = 0.05. |
cc |
logical indicating whether to apply continuity correction, default = FALSE. |
toler |
the desired accuracy using optimise, default = 0.0001. |
logplot |
plot vertical axis as log likelihood, default = FALSE |
supplot |
set minimum likelihood display value in plot, default = -10 |
verb |
show output, default = TRUE. |
$S.val - support for observed OR from expected.
$df - degrees of freedom.
$exp.OR - expected OR.
$S.exp.ORvsObs - support for expected OR versus observed.
$S.exp.ORvsNull - support for expected OR versus the null.
$HAc - Haldane-Anscombe correction applied when a count is 0.
$L.int - likelihood interval of observed OR for specified level of support.
$S_int - specified likelihood interval in units of support.
$observed - observed frequencies.
$expected - the expected values for null hypothesis of no interaction.
$chi.sq - chi-squared statistic.
$corrected - whether chi-squared was corrected, default = FALSE.
$p.value - p value.
$LR.test = the likelihood ratio test statistic.
$lrt.p = the p value for the likelihood ratio test statistic
$residuals - the Pearson residuals.
$alpha - specified significance level.
$conf.int - likelihood-based confidence interval for observed RR.
$all.err.acc - error accuracy for each application of the optimize function.
Aitkin, M. et al (1989) Statistical Modelling in GLIM, Clarendon Press, ISBN : 978-0198522041
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Royall, R. M. (1997). Statistical evidence: A likelihood paradigm. London: Chapman & Hall, ISBN : 978-0412044113
Edwards, A.W.F. (1992) Likelihood, Johns Hopkins Press, ISBN : 978-0801844430
Dienes, Z. (2008) Understanding Psychology as a Science: An Introduction to Scientific and Statistical Inference, Palgrave, MacMillan, ISBN : 978-0230542303
# for folic acid and neural tube defects example, p 146 tab <- as.table(rbind(c(6,587),c(21,581))) dimnames(tab) <- list(Treatment=c("Folic acid","None"),Defect=c("Yes","No")) L_OR(tab, exp.OR = 0.5, L.int = 2, alpha=0.05, cc=FALSE, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)# for folic acid and neural tube defects example, p 146 tab <- as.table(rbind(c(6,587),c(21,581))) dimnames(tab) <- list(Treatment=c("Folic acid","None"),Defect=c("Yes","No")) L_OR(tab, exp.OR = 0.5, L.int = 2, alpha=0.05, cc=FALSE, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)
This function calculates the supports for different regression fits from 2 vectors of data. Models include linear, quadratic and cubic (given sufficient data). A plot is included showing linear (black), quadratic (red) and cubic (blue dashed) lines. P values for the model fits are also given.
L_regress(y, x, verb=TRUE)L_regress(y, x, verb=TRUE)
y |
a numeric vector the same length as x. |
x |
a numeric vector. |
verb |
show output, default = TRUE. |
$S.LNc - corrected support for linear versus null model.
$S.LN - uncorrected support for linear versus null model.
$S.QLc - corrected support for quadratic versus linear model.
$S.QL - uncorrected support for quadratic versus linear model.
S.QCc = support for quadratic versus cubic model.
$N - sample size.
$p.vals - p values for 3 fits.
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Royall, R. M. (1997). Statistical evidence: A likelihood paradigm. London: Chapman & Hall, ISBN : 978-0412044113
Edwards, A.W.F. (1992) Likelihood, Johns Hopkins Press, ISBN : 978-0801844430
# for women's world record times for 1500m event example, p 108 years <- c(0.0, 7.1, 8.9, 8.9, 10.1, 12.8, 17.0, 19.1, 25.0, 28.7, 29.7, 29.9, 35.3, 39.8, 40.2, 41.9, 42.1, 44.0, 44.9, 45.0, 45.1, 45.1, 48.9, 52.9, 53.0, 66.1, 87.9) time <- c(5.30, 5.12, 5.03, 4.79, 4.75, 4.70, 4.63, 4.63, 4.62, 4.59, 4.50, 4.50, 4.32, 4.29, 4.26, 4.21, 4.18, 4.16, 4.12, 4.11, 4.09, 4.02, 3.93, 3.92, 3.87, 3.84, 3.83) m=L_regress(time, years) m# for women's world record times for 1500m event example, p 108 years <- c(0.0, 7.1, 8.9, 8.9, 10.1, 12.8, 17.0, 19.1, 25.0, 28.7, 29.7, 29.9, 35.3, 39.8, 40.2, 41.9, 42.1, 44.0, 44.9, 45.0, 45.1, 45.1, 48.9, 52.9, 53.0, 66.1, 87.9) time <- c(5.30, 5.12, 5.03, 4.79, 4.75, 4.70, 4.63, 4.63, 4.62, 4.59, 4.50, 4.50, 4.32, 4.29, 4.26, 4.21, 4.18, 4.16, 4.12, 4.11, 4.09, 4.02, 3.93, 3.92, 3.87, 3.84, 3.83) m=L_regress(time, years) m
This function calculates the support for an RR from a 2 x 2 categorical data table. An expected RR can be specified and the support calculated for this relative to the observed and null (which is assumed to be 1, but can also be specified) values. A likelihood function is plotted for the obtained RR with a likelihood interval, and expected RR, if specified. The log likelihood plot can optionally be given instead. Chi-squared statistics are also provided and a likelihood-based % confidence interval. It uses the optimize function to locate desired limits for both intervals and other support calculations.
L_RR(table, null=1, exp.RR=NULL, L.int=2, alpha=0.05, cc=FALSE, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)L_RR(table, null=1, exp.RR=NULL, L.int=2, alpha=0.05, cc=FALSE, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)
table |
a 2 x 2 matrix or contingency table containing counts. |
null |
the value against which the obtained RR is tested, default = 1. |
exp.RR |
an expected or hypothetical RR. |
L.int |
likelihood interval given as support values, e.g. 2 or 3, default = 2. |
alpha |
the significance level used, 1 - alpha interval calculated, default = 0.05. |
cc |
logical indicating whether to apply continuity correction, default = FALSE. |
toler |
the desired accuracy using optimise, default = 0.0001. |
logplot |
plot vertical axis as log likelihood, default = FALSE |
supplot |
set minimum likelihood display value in plot, default = -10 |
verb |
show output, default = TRUE. |
$S.val - support for observed RR from expected.
$df - degrees of freedom.
$exp.RR - expected RR.
$S.exp.RRvsObs - support for expected RR versus observed.
$S.exp.RRvsNull - support for expected RR versus the null.
$L.int - likelihood interval of observed RR for specified level of support.
$S_int - specified likelihood interval in units of support.
$observed - observed frequencies.
$expected - the expected values for null hypothesis of no interaction.
$chi.sq - chi-squared statistic.
$corrected - whether chi-squared was corrected, default = FALSE.
$p.value - p value.
$residuals - the Pearson residuals.
$alpha - specified significance level.
$conf.int - likelihood-based confidence interval for observed RR.
$all.err.acc - error accuracy for each application of the optimize function.
Aitkin, M. et al (1989) Statistical Modelling in GLIM, Clarendon Press, ISBN : 978-0198522041
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Royall, R. M. (1997). Statistical evidence: A likelihood paradigm. London: Chapman & Hall, ISBN : 978-0412044113
Edwards, A.W.F. (1992) Likelihood, Johns Hopkins Press, ISBN : 978-0801844430
Dienes, Z. (2008) Understanding Psychology as a Science: An Introduction to Scientific and Statistical Inference, Palgrave, MacMillan, ISBN : 978-0230542303
# for folic acid and neural tube defects example tab <- as.table(rbind(c(6,587),c(21,581))) dimnames(tab) <- list(Treatment=c("Folic acid","None"),Defect=c("Yes","No")) L_RR(tab, exp.RR = 0.5, L.int = 2, alpha=0.05, cc=FALSE, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE) # S. Korea COVID-19 patient mortality tabcor <- as.table(rbind(c(41,3095),c(34,4992))) dimnames(tabcor) <- list(Sex=c("Male","Female"),Status=c("Dead","Alive")) L_RR(tabcor, exp.RR = 0.5, L.int = 2, alpha=0.05, cc=FALSE, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)# for folic acid and neural tube defects example tab <- as.table(rbind(c(6,587),c(21,581))) dimnames(tab) <- list(Treatment=c("Folic acid","None"),Defect=c("Yes","No")) L_RR(tab, exp.RR = 0.5, L.int = 2, alpha=0.05, cc=FALSE, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE) # S. Korea COVID-19 patient mortality tabcor <- as.table(rbind(c(41,3095),c(34,4992))) dimnames(tabcor) <- list(Sex=c("Male","Female"),Status=c("Dead","Alive")) L_RR(tabcor, exp.RR = 0.5, L.int = 2, alpha=0.05, cc=FALSE, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)
This function calculates the required sample size for t tests. The standard deviation and effect size are specified. Calculations given for one sample and independent samples t tests. For a related samples test calculation use the sd for paired differences.
L_t_test_sample_size(MW = 0.05, sd = 1, d = 1.2, S = 3, paired = FALSE, verb=TRUE)L_t_test_sample_size(MW = 0.05, sd = 1, d = 1.2, S = 3, paired = FALSE, verb=TRUE)
MW |
set M1 + W1 probability, default = .05. |
sd |
set standard deviation, default = 1. |
d |
set desired effect size, default = 1.2. |
S |
set strength of evidence (support), default = 3. |
paired |
set to TRUE for one sample and FALSE for independent samples, default = FALSE. |
verb |
show output, default = TRUE. |
$N - required sample size.
$S - specified strength (support) for evidence from the test.
$sd - specified standard deviation.
$d - Cohen's effect size specified.
$m1.w1 - specified probability for combined misleading and weak evidence.
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Cahusac, P.M.B. & Mansour, S.E. (2022) Estimating sample sizes for evidential t tests, Research in Mathematics, 9(1):1-12 https://doi.org/10.1080/27684830.2022.2089373
Royall, R. (2000). "On the Probability of Observing Misleading Statistical Evidence." Journal of the American Statistical Association 95(451): 760.
Royall, R. (2004). The Likelihood paradigm for statistical evidence. The Nature of Scientific Evidence. M. L. Taper and S. R. Lele. Chicago, University of Chicago: 119.
Royall, R. M. (1997). Statistical evidence: A likelihood paradigm. London: Chapman & Hall, ISBN : 978-0412044113
Edwards, A.W.F. (1992) Likelihood, Johns Hopkins Press, ISBN : 978-0801844430
# for one sample or related samples (differences) v = L_t_test_sample_size(MW = 0.2, sd = 1, d = 1, S = 3, paired = TRUE) v # for 2 independent samples v = L_t_test_sample_size(MW = 0.05, sd = 1, d = 1.2, S = 3, paired = FALSE) v# for one sample or related samples (differences) v = L_t_test_sample_size(MW = 0.2, sd = 1, d = 1, S = 3, paired = TRUE) v # for 2 independent samples v = L_t_test_sample_size(MW = 0.05, sd = 1, d = 1.2, S = 3, paired = FALSE) v
This function calculates several different supports. Effect size (Cohen's d) and a second alternative hypothesis value can be specified. The maximum support is the support for the observed mean versus the null value. The support for the specified d versus the null is also calculated. If a second hypothesis value is specified (in units of the original measurements) then two further supports are calculated: d versus 2nd alternative hypothesis, and 2nd alternative hypothesis versus the null. The likelihood curve graphic with MLE and specified hypothesis values is produced. The requested likelihood interval is provided and displayed on likelihood curve. The t and p values for the test against the null value are given.
L_ttest(data1, data2, null=0, d=0.5, alt.2=NULL, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)L_ttest(data1, data2, null=0, d=0.5, alt.2=NULL, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)
data1 |
a (non-empty) numeric vector of data values. |
data2 |
a (non-empty) numeric vector of data values for related sample, default = NULL. |
null |
value for the null hypothesis, default = 0. |
d |
Cohen's effect size, default = 0.5. |
alt.2 |
value for an alternative hypothesis, in units used for data, default = NULL. |
L.int |
likelihood interval given for a given support value, e.g. 2 or 3, default = 2. |
toler |
the desired accuracy using optimise, default = 0.0001. |
logplot |
plot vertical axis as log likelihood, default = FALSE |
supplot |
set minimum likelihood display value in plot, default = -10 |
verb |
show output, default = TRUE. |
$obs.mean - the observed mean or difference in mean for related samples.
$df - degrees of freedom.
$alt.H1 - mean value according to specified d.
$alt.H2 - specified second hypothesis value.
$S_max - maximum support for observed mean against the null.
$S_10 - support for d versus null.
$S_12 - support for d versus specified second hypothesis.
$S_20 - support for second hypothesis versus the null.
$like.int - likelihood interval.
$L.int.spec - specified likelihood interval in units of support.
$null.value - null value.
$t.val - t value for test against null.
$p.val - p value for test against null.
$d.obs - observed effect size.
Cahusac, P.M.B. (2020) Evidence-Based Statistics, Wiley, ISBN : 978-1119549802
Baguley, T. (2012) Serious Stats, Palgrave Macmillan, ISBN: 978-0230577183
Edwards, A.W.F. (1992) Likelihood, Johns Hopkins Press, ISBN : 978-0801844430
Royall, R. M. (2000). On the probability of observing misleading statistical evidence. Journal of the American Statistical Association, 95, 760.
# one sample Gosset's original additional hours of sleep data, p 29 mysample <- c(0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0) L_ttest(mysample, d=.5, alt.2=2, L.int=2) # related samples, p 56 mysample2 <- c(1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4) L_ttest(mysample, mysample2, d=1, alt.2=2, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)# one sample Gosset's original additional hours of sleep data, p 29 mysample <- c(0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0) L_ttest(mysample, d=.5, alt.2=2, L.int=2) # related samples, p 56 mysample2 <- c(1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4) L_ttest(mysample, mysample2, d=1, alt.2=2, L.int=2, toler=0.0001, logplot=FALSE, supplot=-10, verb=TRUE)