> attach (TitanicSurvival) > summary (TitanicSurvival) survived sex age passengerClass no :809 female:466 Min. : 0.1667 1st:323 yes:500 male :843 1st Qu.:21.0000 2nd:277 Median :28.0000 3rd:709 Mean :29.8811 3rd Qu.:39.0000 Max. :80.0000 NA's :263 > ### GOODNESS OF FIT > sextable <- with (TitanicSurvival, table(sex)) > sextable sex female male 466 843 > exp <- c(327, 982) > chisq.test (sextable, p = exp, rescale.p = TRUE) Chi-squared test for given probabilities data: sextable X-squared = 78.761, df = 1, p-value < 2.2e-16 > chisq.test (sextable) Chi-squared test for given probabilities data: sextable X-squared = 108.58, df = 1, p-value < 2.2e-16 > #### BINOMIAL TEST > binom.test(3, 20, p = 0.5) Exact binomial test data: 3 and 20 number of successes = 3, number of trials = 20, p-value = 0.002577 alternative hypothesis: true probability of success is not equal to 0.5 95 percent confidence interval: 0.03207094 0.37892683 sample estimates: probability of success 0.15 ####### FREQUENCIES CALCULATION > titanictable <- table (sex, survived) > titanictable survived sex no yes female 127 339 male 682 161 #### chi-square test for two-way table > chisq.test (titanictable) Pearson's Chi-squared test with Yates' continuity correction data: titanictable X-squared = 363.62, df = 1, p-value < 2.2e-16 > titanictable2 <- table (sex, passengerClass) > titanictable2 passengerClass sex 1st 2nd 3rd female 144 106 216 male 179 171 493 > chisq.test (titanictable2) Pearson's Chi-squared test data: titanictable2 X-squared = 20.379, df = 2, p-value = 0.00003757 > library(gmodels, pos=17) > CrossTable (survived, passengerClass) Cell Contents |-------------------------| | N | | Chi-square contribution | | N / Row Total | | N / Col Total | | N / Table Total | |-------------------------| Total Observations in Table: 1309 | passengerClass survived | 1st | 2nd | 3rd | Row Total | -------------|-----------|-----------|-----------|-----------| no | 123 | 158 | 528 | 809 | | 29.411 | 1.017 | 18.411 | | | 0.152 | 0.195 | 0.653 | 0.618 | | 0.381 | 0.570 | 0.745 | | | 0.094 | 0.121 | 0.403 | | -------------|-----------|-----------|-----------|-----------| yes | 200 | 119 | 181 | 500 | | 47.587 | 1.645 | 29.788 | | | 0.400 | 0.238 | 0.362 | 0.382 | | 0.619 | 0.430 | 0.255 | | | 0.153 | 0.091 | 0.138 | | -------------|-----------|-----------|-----------|-----------| Column Total | 323 | 277 | 709 | 1309 | | 0.247 | 0.212 | 0.542 | | -------------|-----------|-----------|-----------|-----------| ### 3-way table > titanictable3 <- table (sex, survived, passengerClass) > titanictable3 , , passengerClass = 1st survived sex no yes female 5 139 male 118 61 , , passengerClass = 2nd survived sex no yes female 12 94 male 146 25 , , passengerClass = 3rd survived sex no yes female 110 106 male 418 75 > ##### Fishers exact test for 2x2 tables > fisher.test (titanictable) Fisher's Exact Test for Count Data data: titanictable p-value < 2.2e-16 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.06712075 0.11644002 sample estimates: odds ratio 0.08867006 > ##### Fishers exact test for axb tables > fisher.test (titanictable2) Fisher's Exact Test for Count Data data: titanictable2 p-value = 0.00003941 alternative hypothesis: two.sided > ##### Cochran-Mantel-Haenszel a x b x c test > mantelhaen.test (titanictable3) Mantel-Haenszel chi-squared test with continuity correction data: titanictable3 Mantel-Haenszel X-squared = 347.66, df = 1, p-value < 2.2e-16 alternative hypothesis: true common odds ratio is not equal to 1 95 percent confidence interval: 0.06305157 0.11168953 sample estimates: common odds ratio 0.08391782 > cochran <- readXL("D:/Nina/statistica/stat 2019-2020/занятие 7/frq.xlsx", rownames=FALSE, header=TRUE, na="", sheet="Лист2", stringsAsFactors=TRUE) > summary (cochran) pupil result exam_ID a : 3 Min. :0.0 e1:10 b : 3 1st Qu.:0.0 e2:10 c : 3 Median :1.0 e3:10 d : 3 Mean :0.7 e : 3 3rd Qu.:1.0 f : 3 Max. :1.0 (Other):12 > cochran.qtest(result ~ exam_ID|pupil, cochran) Cochran's Q test data: result by exam_ID, block = pupil Q = 5.25, df = 2, p-value = 0.07244 alternative hypothesis: true difference in probabilities is not equal to 0 sample estimates: proba in group e1 proba in group e2 proba in group e3 0.4 0.8 0.9 > ####### LOGISTIC REGRESSION > survival.model <- glm (survived ~ age, family = binomial(), data = TitanicSurvival) > summary (survival.model) Call: glm(formula = survived ~ age, family = binomial(), data = TitanicSurvival) Deviance Residuals: Min 1Q Median 3Q Max -1.1189 -1.0361 -0.9768 1.3187 1.5162 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.136531 0.144715 -0.943 0.3455 age -0.007899 0.004407 -1.792 0.0731 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1414.6 on 1045 degrees of freedom Residual deviance: 1411.4 on 1044 degrees of freedom (263 observations deleted due to missingness) AIC: 1415.4 Number of Fisher Scoring iterations: 4 ########################## SCRIPT############################# attach (TitanicSurvival) summary (TitanicSurvival) ### GOODNESS OF FIT sextable <- with (TitanicSurvival, table(sex)) sextable exp <- c(327, 982) chisq.test (sextable, p = exp, rescale.p = TRUE) chisq.test (sextable) #### BINOMIAL TEST binom.test(3, 20, p = 0.5) ####### FREQUENCIES CALCULATION titanictable <- table (sex, survived) titanictable #### chi-square test for two-way table chisq.test (titanictable) titanictable2 <- table (sex, passengerClass) titanictable2 chisq.test (titanictable2) library(gmodels, pos=17) CrossTable (survived, passengerClass) ### 3-way table titanictable3 <- table (sex, survived, passengerClass) titanictable3 ##### Fishers exact test for 2x2 tables fisher.test (titanictable) ##### Fishers exact test for axb tables fisher.test (titanictable2) ##### Cochran-Mantel-Haenszel a x b x c test mantelhaen.test (titanictable3) cochran <- readXL("D:/Nina/statistica/stat 2019-2020/занятие 7/frq.xlsx", rownames=FALSE, header=TRUE, na="", sheet="Лист2", stringsAsFactors=TRUE) summary (cochran) cochran.qtest(result ~ exam_ID|pupil, cochran) ####### LOGISTIC REGRESSION survival.model <- glm (survived ~ age, family = binomial(), data = TitanicSurvival) summary (survival.model)