Project_3
Siberia
15 марта 2019 г
Intro
Hi! I have chosen as my topic Climate changes and we will explore data about Switzerland. What is behind this decision? Well, I’m fond of a fancy field, namely STS. And Climate changes is one of the disputable topic in the world and it’s one of the reason why STS researchers are into it. I know nothing about Climate changes and decide that with this project I may find the way into it. And I have chosen Switzerland because the UN Environment Regional office for Europe is located in Geneva, but I’m far not sure that I will stick to Switzerland in future.
For today project I will mostly discuss measure inprdsc. The question for this variable from 2012 to 2016
How many people, if any, are there with whom you can discuss intimate and personal matters?
What is behind this decision? Well, recently I’ve read an article called Social Isolation in America: Changes in Core Discussion Networks over Two Decades authors used GSS data to describe how the discussion with close confidants have changed over two decades. They had found that people tend to have less confidants and the proportion of those who had none changed from 10% in 1985 to 24.6% in 2004(well they had more complicated results and more precise questions). So, in my previous project in shadow I’ve take how they answer this question in ESS data across countries.
Intro_intro
As you might see, the proportion for each answer are different across countries, some countries tend to have more None answers than expected, I’ve decide to take a look at one country and find what could explain differences in answers on the question. I’ve chosen Czech and my first idea was to compare these groups in the number of people in household. But it’s count data with low range and not normally distributed, I was trying to find the solution with this variable, but Google highly suggested me to use Poisson distribution which is not a part of our course.
I found more or less continuous variable, which is nice! But it’s highly skewed, but what is important here I haven’t deleted those guys that hadn’t answer to this question. It’s hard to see here, but None group and 1 group are more than other groups haven’t answer to this question(I haven’t conduct test, but it helps we to go further).
On the next step, I’ve decided to compare ages across my groups and this graph starts to look more intuitive and nice despite the fact that bars cross.
Assumptions first step or Pain №2
Okay, now we might be ready to go further but let’s check assumptions first
describeBy(project3_data_ex$agea, project3_data_ex$inprdsc)
##
## Descriptive statistics by group
## group: None
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 242 53.62 14.77 58.5 54.81 12.6 15 92 77 -0.64 -0.22
## se
## X1 0.95
## --------------------------------------------------------
## group: 1
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 527 49.34 16.56 49 49.5 20.76 15 91 76 -0.04 -0.86
## se
## X1 0.72
## --------------------------------------------------------
## group: 2
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 606 45.36 16.26 44 45.3 20.76 15 90 75 0.07 -0.93
## se
## X1 0.66
## --------------------------------------------------------
## group: 3
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 461 44.43 17.76 44 44.21 22.24 15 86 71 0.1 -1.06
## se
## X1 0.83
## --------------------------------------------------------
## group: 4-6
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 299 39.8 17.24 37 38.76 20.76 15 83 68 0.44 -0.85 1
## --------------------------------------------------------
## group: 7-9
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 35 41.03 15.3 41 40.28 17.79 16 73 57 0.43 -0.71
## se
## X1 2.59
## --------------------------------------------------------
## group: 10 or more
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 32 34.28 16.83 30.5 32.62 13.34 15 76 61 0.8 -0.47
## se
## X1 2.97
with(project3_data_ex, by(agea, inprdsc, skewness))
## inprdsc: None
## [1] -0.647571
## --------------------------------------------------------
## inprdsc: 1
## [1] -0.03944026
## --------------------------------------------------------
## inprdsc: 2
## [1] 0.07506741
## --------------------------------------------------------
## inprdsc: 3
## [1] 0.1021307
## --------------------------------------------------------
## inprdsc: 4-6
## [1] 0.4399216
## --------------------------------------------------------
## inprdsc: 7-9
## [1] 0.4474108
## --------------------------------------------------------
## inprdsc: 10 or more
## [1] 0.8340344
with(project3_data_ex, by(agea, inprdsc, kurtosis))
## inprdsc: None
## [1] 2.806537
## --------------------------------------------------------
## inprdsc: 1
## [1] 2.1458
## --------------------------------------------------------
## inprdsc: 2
## [1] 2.078773
## --------------------------------------------------------
## inprdsc: 3
## [1] 1.945645
## --------------------------------------------------------
## inprdsc: 4-6
## [1] 2.164921
## --------------------------------------------------------
## inprdsc: 7-9
## [1] 2.42505
## --------------------------------------------------------
## inprdsc: 10 or more
## [1] 2.696488
#Levene test
leveneTest(project3_data_ex$agea, project3_data_ex$inprdsc, center = median)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 6 4.5233 0.000145 ***
## 2195
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Well, data is not skewed, but different kurtosis functions provide different results and it’s hard to say which results we supposed to use and what range we should apply (-1,1)/(-2,2)/(-3/3). But let’s left it and come back to it later.
But what is not okay, is our Levene test shows us that we violate assumption of homogeneity of variance. It’s not that problem because we can use Welch’s function oneway.test(outcome ~ predictor, data, var.equal = F) which can neglect this problem, but I don’t have enough nerve and I further discover that they have different output(at least in names).
What I’ve decided to do is to choose not that radical country in terms of chi-square results. And this country was Russia.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 6 1.9497 0.06953 .
## 2313
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
It has the same problem with kurtosis, but at least we passed Levene test and we have a little bit better sample sizes. But, I will keep doing the same tests for Czechia, just for the sake of where it could lead us and don’t stop on its results.
ANOVA
Our hypothesis:
- H0: There is no difference in means in all groups.
- H1: There is significant difference in means at least in one group.
Ru anova
anova_model <- aov(agea ~ inprdsc, project3_data_2)
summary(anova_model)
## Df Sum Sq Mean Sq F value Pr(>F)
## inprdsc 6 18038 3006.3 9.402 3.22e-10 ***
## Residuals 2313 739568 319.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Our ratio of systematic variance to unsystematic variance or simply F-ratio is 9.402 which in terms of p.value is significant.
Ru come back to assumptions
plotNormalHistogram(residuals(anova_model))
On the first graph we could see the variance in each group,as long as this graph don’t take funnel shape we are fine(I guess) and variance across groups are similar which was also the conclusion of Levene test. From the second graph we could estimate the normality of residuals in the model. In our case it doesn’t look normal. On the third graph we can see distribution of our residuals and compare with blue line. Well, it’s not good but I guess it could be much worse. Nevertheless, we have to check our results with non-parametric test or something more robust.
Cz Welch’s Anova
##
## One-way analysis of means (not assuming equal variances)
##
## data: agea and inprdsc
## F = 24.19, num df = 6.00, denom df = 232.73, p-value < 2.2e-16
For some reason aov and oneway.test have different output and I cannot plot these functions plot(anova_model_ex) plotNormalHistogram(residuals(anova_model_ex))
Ru post hoc Bonferroni and Benjamini–Hochberg
Let’s first assume that we have no problem with normality and we can verify our results with post hoc tests.
pairwise.t.test(project3_data_2$agea, project3_data_2$inprdsc,
p.adjust.method = "bonferroni", paired = F)
##
## Pairwise comparisons using t tests with pooled SD
##
## data: project3_data_2$agea and project3_data_2$inprdsc
##
## None 1 2 3 4-6 7-9
## 1 0.04209 - - - - -
## 2 7.3e-07 0.01439 - - - -
## 3 6.6e-09 0.00013 1.00000 - - -
## 4-6 2.4e-06 0.02078 1.00000 1.00000 - -
## 7-9 0.08102 1.00000 1.00000 1.00000 1.00000 -
## 10 or more 0.10478 1.00000 1.00000 1.00000 1.00000 1.00000
##
## P value adjustment method: bonferroni
pairwise.t.test(project3_data_2$agea, project3_data_2$inprdsc,
p.adjust.method = "BH", paired = F)
##
## Pairwise comparisons using t tests with pooled SD
##
## data: project3_data_2$agea and project3_data_2$inprdsc
##
## None 1 2 3 4-6 7-9
## 1 0.0060 - - - - -
## 2 3.6e-07 0.0029 - - - -
## 3 6.6e-09 3.4e-05 0.2664 - - -
## 4-6 8.0e-07 0.0035 0.7134 0.8592 - -
## 7-9 0.0101 0.2664 0.8822 0.9035 0.9911 -
## 10 or more 0.0116 0.3764 0.9617 0.7134 0.8592 0.8822
##
## P value adjustment method: BH
We use post hoc test because we don’t have any specific hypothesis, we compare each group with other and trying to reveal which group means differ, it’s like t-test but we don’t make family-wise error by multiple comparison(as long as we don’t violate assumptions especially equivalence of group variances).
First results are from Bonferroni’s test which controls the Type 1 error rate and it’s kinda conservative in terms of we have a trade-off with higher chance of Type 2 error. Simply speaking we divide Type 1 error rate by the number of comparison and it’s our critical value for each pair of groups. So, the more group we have the smaller chance we have to reveal any difference between groups. But we can be more or less sure in our results. In my case there are 7 pairs of groups which means are different.
Second results are from Benjamini-Hochberg’s test which also controls the Type 1 error but It’s less conservative and take more step-up way of measuring critical value. In my case there are 9 pairs of groups which means are different.
Ru post hoc Tukey
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = agea ~ inprdsc, data = project3_data_2)
##
## $inprdsc
## diff lwr upr p adj
## 1-None -4.56397083 -8.918177 -0.2097650 0.0328257
## 2-None -7.98590041 -12.243637 -3.7281642 0.0000007
## 3-None -9.54973361 -14.009142 -5.0903251 0.0000000
## 4-6-None -8.91190073 -13.856179 -3.9676221 0.0000024
## 7-9-None -8.87829640 -17.936101 0.1795078 0.0590094
## 10 or more-None -7.71930814 -15.824557 0.3859408 0.0738945
## 2-1 -3.42192959 -6.391816 -0.4520427 0.0121569
## 3-1 -4.98576278 -8.238183 -1.7333430 0.0001298
## 4-6-1 -4.34792991 -8.238663 -0.4571971 0.0171463
## 7-9-1 -4.31432557 -12.842850 4.2141988 0.7492983
## 10 or more-1 -3.15533731 -10.664463 4.3537886 0.8784767
## 3-2 -1.56383319 -4.685923 1.5582565 0.7580042
## 4-6-2 -0.92600032 -4.708461 2.8564607 0.9912604
## 7-9-2 -0.89239598 -9.372074 7.5872819 0.9999279
## 10 or more-2 0.26659227 -7.187010 7.7201942 0.9999999
## 4-6-3 0.63783287 -3.370286 4.6459518 0.9992033
## 7-9-3 0.67143721 -7.911275 9.2541491 0.9999875
## 10 or more-3 1.83042547 -5.740188 9.4010389 0.9918253
## 7-9-4-6 0.03360434 -8.810745 8.8779534 1.0000000
## 10 or more-4-6 1.19259259 -6.673395 9.0585798 0.9993964
## 10 or more-7-9 1.15898826 -9.771322 12.0892984 0.9999246
postHocs <- glht(anova_model, linfct = mcp(inprdsc = "Tukey"))
summary(postHocs)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = agea ~ inprdsc, data = project3_data_2)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 - None == 0 -4.5640 1.4755 -3.093 0.0276 *
## 2 - None == 0 -7.9859 1.4428 -5.535 <0.001 ***
## 3 - None == 0 -9.5497 1.5112 -6.319 <0.001 ***
## 4-6 - None == 0 -8.9119 1.6755 -5.319 <0.001 ***
## 7-9 - None == 0 -8.8783 3.0695 -2.892 0.0498 *
## 10 or more - None == 0 -7.7193 2.7467 -2.810 0.0622 .
## 2 - 1 == 0 -3.4219 1.0064 -3.400 0.0104 *
## 3 - 1 == 0 -4.9858 1.1022 -4.524 <0.001 ***
## 4-6 - 1 == 0 -4.3479 1.3185 -3.298 0.0143 *
## 7-9 - 1 == 0 -4.3143 2.8901 -1.493 0.7206
## 10 or more - 1 == 0 -3.1553 2.5447 -1.240 0.8610
## 3 - 2 == 0 -1.5638 1.0580 -1.478 0.7299
## 4-6 - 2 == 0 -0.9260 1.2818 -0.722 0.9895
## 7-9 - 2 == 0 -0.8924 2.8736 -0.311 0.9999
## 10 or more - 2 == 0 0.2666 2.5259 0.106 1.0000
## 4-6 - 3 == 0 0.6378 1.3583 0.470 0.9990
## 7-9 - 3 == 0 0.6714 2.9085 0.231 1.0000
## 10 or more - 3 == 0 1.8304 2.5655 0.713 0.9902
## 7-9 - 4-6 == 0 0.0336 2.9971 0.011 1.0000
## 10 or more - 4-6 == 0 1.1926 2.6656 0.447 0.9993
## 10 or more - 7-9 == 0 1.1590 3.7040 0.313 0.9999
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
plot(TukeyHSD(anova_model))
Another popular post hoc test is Tukey test, I don’t know why it’s that good or popular the only thing I remember it doesn’t work well with different group sizes and unequal variances. I performed 2 different functions one is from already installed in R package and another from gmodels package. They give slightly different results and I don’t know why. But let’s work with TukeyHSD. The results in terms of pairs of groups are the same as Bonferroni. On the graph we can see all pairs and their confidence intervals and which of them don’t cross zero point, they are signal to us about differences in means. We could also say that people who have no one to speak about personal matters are usually elder than people who have from 1 to 6 close confidants. And it’s also works for people who have only 1 person.
Cz post hoc test
##
## Pairwise comparisons using t tests with pooled SD
##
## data: project3_data_ex$agea and project3_data_ex$inprdsc
##
## None 1 2 3 4-6 7-9
## 1 0.01963 - - - - -
## 2 1.7e-09 0.00130 - - - -
## 3 9.5e-11 8.3e-05 1.00000 - - -
## 4-6 < 2e-16 7.7e-14 4.9e-05 0.00378 - -
## 7-9 0.00062 0.08937 1.00000 1.00000 1.00000 -
## 10 or more 1.6e-08 1.5e-05 0.00517 0.01803 1.00000 1.00000
##
## P value adjustment method: bonferroni
##
## Pairwise comparisons using t tests with pooled SD
##
## data: project3_data_ex$agea and project3_data_ex$inprdsc
##
## None 1 2 3 4-6 7-9
## 1 0.00140 - - - - -
## 2 4.3e-10 0.00013 - - - -
## 3 3.2e-11 1.0e-05 0.38382 - - -
## 4-6 < 2e-16 3.8e-14 7.0e-06 0.00034 - -
## 7-9 6.9e-05 0.00596 0.15665 0.26926 0.67830 -
## 10 or more 3.2e-09 2.5e-06 0.00043 0.00139 0.09813 0.12021
##
## P value adjustment method: BH
I cannot produce Tukey test because of wrong input format
Ru non-parametric + Dunn post hoc
kruskal.test(agea ~ inprdsc, project3_data_2)
##
## Kruskal-Wallis rank sum test
##
## data: agea by inprdsc
## Kruskal-Wallis chi-squared = 51.309, df = 6, p-value = 2.567e-09
DunnTest(agea ~ inprdsc, project3_data_2, method = "bonferroni")
##
## Dunn's test of multiple comparisons using rank sums : bonferroni
##
## mean.rank.diff pval
## 1-None -165.227406 0.05867 .
## 2-None -287.763107 2.1e-06 ***
## 3-None -341.506466 3.4e-08 ***
## 4-6-None -318.336623 8.2e-06 ***
## 7-9-None -323.850625 0.10186
## 10 or more-None -275.919957 0.15371
## 2-1 -122.535701 0.02419 *
## 3-1 -176.279060 0.00041 ***
## 4-6-1 -153.109218 0.04060 *
## 7-9-1 -158.623219 1.00000
## 10 or more-1 -110.692551 1.00000
## 3-2 -53.743358 1.00000
## 4-6-2 -30.573516 1.00000
## 7-9-2 -36.087518 1.00000
## 10 or more-2 11.843151 1.00000
## 4-6-3 23.169842 1.00000
## 7-9-3 17.655841 1.00000
## 10 or more-3 65.586509 1.00000
## 7-9-4-6 -5.514002 1.00000
## 10 or more-4-6 42.416667 1.00000
## 10 or more-7-9 47.930668 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Well, we return to try non-parametric tests. We performed Kruskal-Wallis test and we have a significant results, which means that it’s greater than critical chi-square value and there is difference at least in one group from another in terms of median.
Dunn post hoc test gives us almost the same results as Tukey test, but difference between None and 1 group is not significant now, but we again used Bonferroni method which could be too harsh for our model.
Ru bootstrap or Pain №3
t1waybt(agea ~ inprdsc, data = project3_data_2, tr = 0.2, nboot = 5000)
## Call:
## t1waybt(formula = agea ~ inprdsc, data = project3_data_2, tr = 0.2,
## nboot = 5000)
##
## Effective number of bootstrap samples was 5000.
##
## Test statistic: 7.5257
## p-value: 0
## Variance explained 0.046
## Effect size 0.216
mcppb20(agea ~ inprdsc, data = project3_data_2, tr = 0.2, nboot = 5000)
## Call:
## mcppb20(formula = agea ~ inprdsc, data = project3_data_2, tr = 0.2,
## nboot = 5000)
##
## psihat ci.lower ci.upper p-value
## 2 vs. 1 5.79677 -0.28186 12.14714 0.0028
## 2 vs. 3 9.80103 3.47243 15.35625 0.0000
## 2 vs. 4-6 11.31717 5.37078 17.19112 0.0000
## 2 vs. None 10.64716 3.61070 16.87436 0.0000
## 2 vs. 10 or more 10.13210 -1.00538 21.55025 0.0092
## 2 vs. 7-9 8.22269 -2.96639 19.82353 0.0306
## 1 vs. 3 4.00426 0.04066 7.81062 0.0016
## 1 vs. 4-6 5.52040 0.96666 9.70110 0.0004
## 1 vs. None 4.85039 -0.63073 9.79502 0.0068
## 1 vs. 10 or more 4.33533 -5.98559 15.30605 0.2438
## 1 vs. 7-9 2.42592 -8.19461 13.34803 0.4956
## 3 vs. 4-6 1.51614 -2.49421 5.38046 0.2572
## 3 vs. None 0.84614 -4.27493 5.81757 0.6052
## 3 vs. 10 or more 0.33107 -10.34282 10.95341 0.9160
## 3 vs. 7-9 -1.57834 -11.75481 9.20518 0.6784
## 4-6 vs. None -0.67001 -5.79849 4.63001 0.6692
## 4-6 vs. 10 or more -1.18507 -11.45767 9.22055 0.7496
## 4-6 vs. 7-9 -3.09448 -14.06648 7.75604 0.4096
## None vs. 10 or more -0.51506 -11.58667 11.25012 0.8986
## None vs. 7-9 -2.42447 -12.90269 9.16703 0.5408
## 10 or more vs. 7-9 -1.90941 -16.63176 13.27412 0.7108
I also tried to do bootstrap but it’s a mess and I don’t know what’s wrong.
Cz non-parametric + post hoc
##
## Kruskal-Wallis rank sum test
##
## data: agea by inprdsc
## Kruskal-Wallis chi-squared = 132.32, df = 6, p-value < 2.2e-16
##
## Dunn's test of multiple comparisons using rank sums : bonferroni
##
## mean.rank.diff pval
## 1-None -168.20881 0.01376 *
## 2-None -314.34594 1.7e-09 ***
## 3-None -347.85199 1.1e-10 ***
## 4-6-None -522.27093 < 2e-16 ***
## 7-9-None -478.98335 0.00065 ***
## 10 or more-None -723.85925 3.0e-08 ***
## 2-1 -146.13713 0.00239 **
## 3-1 -179.64319 0.00020 ***
## 4-6-1 -354.06212 3.0e-13 ***
## 7-9-1 -310.77455 0.10709
## 10 or more-1 -555.65044 3.3e-05 ***
## 3-2 -33.50605 1.00000
## 4-6-2 -207.92499 7.8e-05 ***
## 7-9-2 -164.63741 1.00000
## 10 or more-2 -409.51330 0.00804 **
## 4-6-3 -174.41894 0.00462 **
## 7-9-3 -131.13136 1.00000
## 10 or more-3 -376.00725 0.02550 *
## 7-9-4-6 43.28758 1.00000
## 10 or more-4-6 -201.58832 1.00000
## 10 or more-7-9 -244.87589 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Ru Effect size
We can also estimate our effect size.
## eta.sq eta.sq.part
## inprdsc 0.02380865 0.02380865
Which is actually quite low, this could mean that despite the fact that we have statistically significant differences in ages among people with different number of confidants with whom they can speak about personal, but the effect of age is far not substantive.
Conclusion
The result could be reported as something like this:
There was a significant effect of age on the number of close confidants for our respondents, F(6,2313) = 9.402, p < 0.05
or
Bonferroni tests revealed significant differences between group who has none confidants and groups with 1-6 confidants, and also group with only one confidant and groups with 2-6 confidants.
To make it may be more clear try not to think about vague term close confidant or personal matter. But as something like strong ties(I’m not sure that it’s good here but whatever). When you growing up this number changes, but after some time you start to lose more ties than gain new. In Czech example we’ve seen that the same pairs of groups are used, but they have even more and it makes harder to explain it intuitively. And it’s also an issue to explain why groups 7-9 and 10 + are not significant or at least not in all tests. May be the problem is a low number of these answers ~50 in each group, but together they represent only 4% of Russian sample.