R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> # data
> 
> x <- 2
> n <- 25
> 
> # mle
> 
> phat <- x / n
> 
> # confidence level
> 
> conf.level <- 0.95
> 
> # Wald interval (used to be standard taught in intro stats)
> 
> phat + c(-1, 1) * qnorm((1 + conf.level) / 2) * sqrt(phat * (1 - phat) / n)
[1] -0.02634498  0.18634498
> 
> # score interval (some intro stats books now teach this)
> 
> prop.test(x, n, conf.level = conf.level, correct = FALSE)

	1-sample proportions test without continuity correction

data:  x out of n, null probability 0.5
X-squared = 17.64, df = 1, p-value = 2.669e-05
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.0222204 0.2496611
sample estimates:
   p 
0.08 

> 
> # the "correct = FALSE" is just bizarre.  (Does that mean we don't want the
> # correct answer?)  No.  The R statement help(prop.test) explains that it
> # it means we do not want to use "continuity correction".  And we don't.
> # No authority recommends what prop.test does by default.
> # Agresti, Section 1.4.2, recommends correct = FALSE.
> # Brown, Cai and DasGupta (Statistical Science, 2005, 20, pp. 375-379)
> # criticize Geyer and Meeden (Statistical Science, 2005, 20, pp. 358-366)
> # for using prop.test with correct = TRUE, providing plots of coverage
> # probability for with correct = FALSE and correct = TRUE to show this.
> # AFAIK no authority defends correct = TRUE, even help(prop.test) cites
> # no such authority.  They do cite
> # Newcombe (Statistics in Medicine, 1998, 17, 857-872) where correct = FALSE
> # is his method 3 and correct = TRUE is (apparently) his method 4
> # Newcombe's simulations do not show the correct = TRUE is better than
> # correct = FALSE and his conclusion does not recommend correct = TRUE
> # over correct = FALSE (but he does not pick a particular method to recommend).
> #
> # for a full derivation of the score interval and a check that correct = FALSE
> # actually calculates it see
> # http://www.stat.umn.edu/geyer/5102/slides/s2.pdf, slides 113-116 and 123
> # (hmmmm.  That does not actually seem to check.  We'll check here.)
> crit <- qnorm((1 + conf.level) / 2)
> foo <- (phat + crit^2 / (2 * n) + c(-1, 1) * crit *
+     sqrt(phat * (1 - phat) / n + crit^2 / (4 * n^2))) / (1 + crit^2 / n)
> foo
[1] 0.0222204 0.2496611
> 
> # just get the interval
> 
> bar <- prop.test(x, n, conf.level = conf.level, correct = FALSE)
> names(bar)
[1] "statistic"   "parameter"   "p.value"     "estimate"    "null.value" 
[6] "conf.int"    "alternative" "method"      "data.name"  
> bar$conf.int
[1] 0.0222204 0.2496611
attr(,"conf.level")
[1] 0.95
> all.equal(bar$conf.int, foo, check.attributes = FALSE)
[1] TRUE
> 
> # level set of log likelihood interval (too esoteric for intro stats)
> 
> logl <- function(p, x, n) ifelse(x == 0, n * log(1 - p),
+     ifelse(x == n, n * log(p), x * log(p) + (n - x) * log(1 - p)))
> crit <- qchisq(conf.level, df = 1)
> fred <- function(p) 2 * (logl(phat, x, n) - logl(p, x, n)) - crit
> tol <- sqrt(.Machine$double.eps)
> if (phat == 0) {
+     low <- 0
+ } else {
+     low <- uniroot(fred, lower = 0, upper = phat, tol = tol)$root
+ }
> if (phat == 1) {
+     hig <- 1
+ } else {
+     hig <- uniroot(fred, lower = phat, upper = 1, tol = tol)$root
+ }
> c(low, hig)
[1] 0.01376568 0.22711456
> 
> # can also do with confint function in MASS package ???
> 
> library(MASS)
> gout <- glm(rbind(c(x, n - x)) ~ 1, family = binomial)
> cout <- confint(gout, level = conf.level)
Waiting for profiling to be done...
> cout
    2.5 %    97.5 % 
-4.272091 -1.225066 
> fred <- 1 / (1 + exp(- cout))
> fred
     2.5 %     97.5 % 
0.01376058 0.22704616 
> all.equal(c(low, hig), fred, check.attributes = FALSE, tolerance = 4e-4)
[1] TRUE
> 
> # intervals compared
> 
> phat + c(-1, 1) * qnorm((1 + conf.level) / 2) * sqrt(phat * (1 - phat) / n)
[1] -0.02634498  0.18634498
> as.vector(bar$conf.int)
[1] 0.0222204 0.2496611
> c(low, hig)
[1] 0.01376568 0.22711456
> 
> # Hypothesis tests
> 
> # data
> 
> x <- 1300
> n <- 2500
> 
> # mle
> 
> phat <- x / n
> 
> # upper tailed test
> # null hypothesis
> 
> p0 <- 1 / 2
> 
> # test statistic and P-value for score test (standard in intro stats)
> 
> tstat <- (x - n * p0) / sqrt(n * p0 * (1 - p0))
> tstat
[1] 2
> pnorm(tstat, lower.tail = FALSE)
[1] 0.02275013
> 
> # P-value for exact test (should be standard in intro stats ?????)
> pbinom(x - 1, n, p0, lower.tail = FALSE)
[1] 0.02384093
> 
> # fuzzy P-value for exact test (Geyer and Meeden, Statistical Science,
> # 2005, vol 20, pp. 358-387) is uniformly distributed on interval
> pbinom(c(x, x - 1), n, p0, lower.tail = FALSE)
[1] 0.02168093 0.02384093
> 
> # test statistic and P-value for likelihood ratio test
> # for one-tailed test have to use signed likelihood ratio statistic
> 
> tstat <- 2 * (logl(phat, x, n) - logl(p0, x, n))
> tstat <- sqrt(tstat)
> tstat <- tstat * sign(phat - p0)
> tstat
[1] 2.000267
> pnorm(tstat, lower.tail = FALSE)
[1] 0.02273573
> 
> # test statistic and P-value for Wald test
> tstat <- (x - n * p0) / sqrt(n * phat * (1 - phat))
> tstat
[1] 2.001602
> pnorm(tstat, lower.tail = FALSE)
[1] 0.02266378
> 
> # two tailed test (same data)
> 
> tstat <- (x - n * p0) / sqrt(n * p0 * (1 - p0))
> tstat
[1] 2
> pnorm(abs(tstat), lower.tail = FALSE) * 2
[1] 0.04550026
> 
> # P-value for exact test (no exact two-tailed test)
> 
> # fuzzy P-value for exact test (Geyer and Meeden, Statistical Science,
> # 2005, vol 20, pp. 358-387) now not simple
> library(ump)
> print(arpv.binom(x, n, p0, plot = FALSE))
$alpha
[1] 0.04336187 0.04336187 0.04768187

$phi
[1] 0.00000e+00 1.73133e-12 1.00000e+00

> 
> # test statistic and P-value for likelihood ratio test
> tstat <- 2 * (logl(phat, x, n) - logl(p0, x, n))
> pchisq(tstat, df = 1, lower.tail = FALSE)
[1] 0.04547146
> 
> # test statistic and P-value for Wald test
> tstat <- (x - n * p0) / sqrt(n * phat * (1 - phat))
> tstat
[1] 2.001602
> pnorm(abs(tstat), lower.tail = FALSE) * 2
[1] 0.04532756
> 
> 
> proc.time()
   user  system elapsed 
  0.160   0.028   0.359