Hi all, I'm working with hospital-based surveillance data. My peers and I are trained (loosely speaking) as frequentists. I'm mindful of the aphorism:
Statistics: A subject which most statisticians find difficult but which many physicians are experts on.
... and want to improve through practice.
Here's a good example.
hospitalizations <- structure(list(timespan = structure(1:4, .Label = c("2002-2004", "2005-2007", "2008-2010", "2011-2013"), class = c("ordered", "factor")), c(65, 66, 52, 44), c(4417, 4361, 4458, 4560)), class = "data.frame", row.names = c(NA, -4L), .Names = c("timespan", "count", "pop"))
Trend or No Trend?
In this particular case, I've been advised to rely on
stats::prop.trend.test(), which exploits the sequencing of the observations, in preference to
> with(hospitalizations, prop.test(count, pop)) 4-sample test for equality of proportions without continuity correction data: count out of pop X-squared = 7.2, df = 3, p-value = 0.07 alternative hypothesis: two.sided sample estimates: prop 1 prop 2 prop 3 prop 4 0.014716 0.015134 0.011664 0.009649
Performs chi-squared test for trend in proportions, i.e., a test asymptotically optimal for local alternatives where the log odds vary in proportion with
score. By default,
scoreis chosen as the group numbers.
> with(hospitalizations, prop.trend.test(count, pop)) Chi-squared Test for Trend in Proportions data: count out of pop , using scores: 1 2 3 4 X-squared = 6.2, df = 1, p-value = 0.01
Doing it Differently
How should I start learning about ways to tackle this question, and related ones, using a different, complementary approach? I'm particularly interested in using this as a jumping-off point into simple Bayesian / MLE techniques.
For example, suppose I'm willing to start with a prior based on the first two observations, and then update my belief/model according to the next two. How would I implement that in R? Any recommendations for libraries, whitepapers, chapters, books, articles, posts, etc.? (Stan sounds fun.)
Or, suppose I want to generate a maximum-likelihood estimate of a slope and intercept. I don't know where to start. (I'm an R user, and like to learn from code-based examples.)
I recognize that this is a wide-open question, and you have other jobs to do. Please allow me to thank you in advance for your time and attention. Comments on this specific example (
prop.test()) are also very welcome.