SOME STATISTICAL CALCULATORS

Chi square test with counts

Đôi khi trong những tình huống phân tích dữ liệu chúng ta cần kiểm định sự khác biệt số liệu trong các bảng bằng Chi square test.

Chẳng hạn trong bảng sau, chúng ta cần biết khác biệt về giới tính giữa các nhóm có ý nghĩa thống kê hay không.

M <- as.table(rbind(c(732, 425, 317), c(484, 251, 318)))
Xsq <- chisq.test(M)  # Prints test summary
Xsq
Xsq$p.value
Xsq$observed   # observed counts (same as M)
Xsq$expected   # expected counts under the null
Xsq$residuals  # Pearson residuals
Xsq$stdres     # standardized residuals

p = 2.319e-06 < 0.001

Một số công thức khác

## Effect of simulating p-values using Monte Carlo simulation
x <- matrix(c(12, 5, 7, 7), ncol = 2)
chisq.test(x)$p.value           # 0.4233
chisq.test(x, simulate.p.value = TRUE, B = 10000)$p.value # around 0.29!

x <- c(A = 20, B = 15, C = 25)
chisq.test(x)
chisq.test(as.table(x))      

x <- c(89,37,30,28,2)
p <- c(40,20,20,15,5)
try(
  chisq.test(x, p = p)                # gives an error
)
chisq.test(x, p = p, rescale.p = TRUE)

p <- c(0.40, 0.20, 0.20, 0.19, 0.01)

# Expected count in category 5
# is 1.86 < 5 ==> chi square approx.
chisq.test(x, p = p)            #               maybe doubtful, but is ok!
chisq.test(x, p = p, simulate.p.value = TRUE)

prop.test

Khác biệt tỉ lệ giữa các nhóm có ý nghĩa không

Khác biệt số đếm so với tỉ lệ p = 0.87 của dân số

prop.test(x=25,n=48, p=0.87, alternative='less')

prop.test(x=25, n=48, p = 0.87,
          alternative = c("two.sided", "less", "greater"),
          conf.level = 0.95, correct = TRUE)

Fisher Exact test

# Fisher's exact test
x <- matrix(c(22, 42, 16, 1, 98, 98, 98, 98), byrow = TRUE, nrow = 2, ncol = 4);

fisher.test(x)

Khoảng tin cậy 95% của 2 tỉ lệ

library(pgirmess)
x<-c(2,10,7,8,7) # eg: number of positive cases
y<-c(56,22,7,20,5)# eg: number of negative cases

CI(cbind(x,y), conf.level=0.95)

Value

A 3 column matrix.

• Column 1: percentage estimate

• Column 2: lower limit of the confidence interval

• column 3: upper limit of the confidence interval

Cách khác:

x <- matrix(c(10, 15, 27, 31, 10, 18), ncol = 2)
CI(x, conf.level= 0.95)

Tính thủ công

# 95CI for proportion differnce
p1 = 0
p2 = 0.012
n1 = 98
n2 = 150
z = 1.96
(p1 - p2)  -  z * sqrt((p1*(1 - p1)/n1 + p2*(1 - p2)/n2))    # LOW 95CI
(p1 - p2)  +  z * sqrt((p1*(1 - p1)/n1 + p2*(1 - p2)/n2))    # HIGH 95CI

confidence intervals for a binomial proportion

library(DescTools)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

binom.test(7, 21,
           0.5,
           alternative="two.sided",
           conf.level=0.95)

    Exact binomial test

data:  7 and 21
number of successes = 7, number of trials = 21, p-value =
0.1892
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.1458769 0.5696755
sample estimates:
probability of success 
             0.3333333 

observed = c(7, 14, 25, 37)
total = sum(observed)

BinomCI(observed, total,
        conf.level = 0.95,
        method = "clopper-pearson")

           est     lwr.ci    upr.ci
x.1 0.08433735 0.03458229 0.1660766
x.2 0.16867470 0.09539963 0.2667839
x.3 0.30120482 0.20531236 0.4117773
x.4 0.44578313 0.33655876 0.5589902

confidence intervals for a binomial proportion

BinomCI(7, 21,
        conf.level = 0.95,
        method = "clopper-pearson")

           est    lwr.ci    upr.ci
[1,] 0.3333333 0.1458769 0.5696755

CI for proportion difference

library(PropCIs)
# p1 = 7/12, p2 = 13/17
diffscoreci(7, 21, 13, 17,
            conf.level=0.95)

data:  

95 percent confidence interval:
 -0.6685985 -0.1103804

confidence intervals for a multinomial proportion

Sex Count Female 10
Male 9 Other 1 No answer 1 ———- —– Total 21

library(DescTools)

observed = c(10, 9, 1, 1)

MultinomCI(observed,
           conf.level=0.95,
           method="sisonglaz")

            est    lwr.ci    upr.ci
[1,] 0.47619048 0.2857143 0.7009460
[2,] 0.42857143 0.2380952 0.6533270
[3,] 0.04761905 0.0000000 0.2723746
[4,] 0.04761905 0.0000000 0.2723746

CIs for numeric variables

library(epiDisplay)
data(Oswego)
.data <- Oswego
attach(.data)
# logical variable
ci(ill)
# numeric variable
ci(age)
# factor
ci(sex=="M")
ci(sex=="F")
detach(.data)

Thank you for reading