Importing the Data

ESS <- read.spss("ESS7SI.sav", use.value.labels = T, to.data.frame = T)   

data1 <- select(ESS, cgtsday, mnactic, fltdpr, alcwknd, gndr, agea, hlpfmhr)
data1 <- na.omit(data1)

str(data1)

## 'data.frame':    73 obs. of  7 variables:
##  $ cgtsday: Factor w/ 22 levels "1","2","3","4",..: 17 4 7 17 6 10 5 17 17 20 ...
##  $ mnactic: Factor w/ 9 levels "Paid work","Education",..: 8 6 8 1 1 1 1 6 1 1 ...
##  $ fltdpr : Factor w/ 4 levels "None or almost none of the time",..: 1 2 3 3 2 2 1 3 1 1 ...
##  $ alcwknd: Factor w/ 87 levels "0","4","7","9",..: 18 5 5 5 15 11 5 8 40 27 ...
##  $ gndr   : Factor w/ 2 levels "Male","Female": 2 2 2 2 2 2 2 2 1 1 ...
##  $ agea   : Factor w/ 78 levels "15","16","17",..: 15 64 34 10 27 43 35 45 25 34 ...
##  $ hlpfmhr: Factor w/ 7 levels "1-10 hours a week",..: 1 1 1 1 2 1 1 1 1 1 ...
##  - attr(*, "variable.labels")= Named chr  "Title of dataset" "ESS round" "Edition" "Production date" ...
##   ..- attr(*, "names")= chr  "name" "essround" "edition" "proddate" ...
##  - attr(*, "codepage")= int 65001
##  - attr(*, "na.action")= 'omit' Named int  1 2 3 4 5 6 7 8 9 10 ...
##   ..- attr(*, "names")= chr  "1" "2" "3" "4" ...

data1$alcwknd <- as.numeric(as.character(data1$alcwknd))
data1$cgtsday <- as.numeric(as.character(data1$cgtsday))
data1$agea <- as.numeric(as.character(data1$agea))

Recoding the Data

Let’s recode the hours helping others variable (hlpfmhr) so the variable observations correspond to the amount of hours.

data1$mnactic <- recode(data1$mnactic, "Housework, looking after children, others" == "Housework")

data1$hlpfmhr <- recode(data1$hlpfmhr, "1" == "1-10")
data1$hlpfmhr <- recode(data1$hlpfmhr, "2" == "11-20")
data1$hlpfmhr <- recode(data1$hlpfmhr, "3" == "21-30")
data1$hlpfmhr <- recode(data1$hlpfmhr, "4" == "31-40")
data1$hlpfmhr <- recode(data1$hlpfmhr, "5" == "41-50")
data1$hlpfmhr <- recode(data1$hlpfmhr, "6" == "55+")
data1$hlpfmhr <- recode(data1$hlpfmhr, "55" == "<1")
data1$hlpfmhr <- recode(data1$hlpfmhr, "hours a week" == "")
data1$hlpfmhr <- recode(data1$hlpfmhr, "hour a week" == "")

Hypotheses

h0 <- as.data.frame(cbind("H0:", "Gender and being exposed to breathing smoke in a job are independent variables."))
h1 <- as.data.frame(cbind("H1:", "Gender and being exposed to breathing smoke in a job are not independent variables."))
Hs <- as.data.frame(rbind(h0, h1))
colnames(Hs) <- c(" ", " ")

Hs %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F) %>%
  row_spec(1:2, background = "White")

Checking Assumptions

sjp.xtab(data2$jbexebs, data2$gndr, bar.pos = "stack", geom.colors = c("skyblue", "salmon"), 
         title = "Being exposed to breathing smoke in a job distributed by gender", axis.titles = "Exposure to breathing smoke in a job", 
         legend.title = "Gender", axis.labels = c("Not exposed", "Exposed"), show.total = FALSE, margin = "row")

From the barplot it can be seen that there are much less people in Slovenia who are exposed to breathing smoke in a job. Nevertheless, there are more men than women among them.

We build contingency table and check the assumptions.

cont_table <- table(data2$gndr, data2$jbexebs, dnn = c("Gender", "Exposed to breathing smoke in a job"))
cont_table

##         Exposed to breathing smoke in a job
## Gender   Not marked Marked
##   Male           92     47
##   Female        134     14

From the contingency table it is seen that we met the assumptions of Chi-squared test:

• There are no empty cells.
  
• Less than 20% of observations are 5 or less.

Test & Residuals

test <- chisq.test(cont_table)  
test

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  cont_table
## X-squared = 20, df = 1, p-value = 1e-06

The X-squared(1) = 70.206 and with p-value < 2.2e-16 we can claim that there is a statistically significant association between breathing smoke in a job and gender.

We analyze the residuals.

test$stdres

##         Exposed to breathing smoke in a job
## Gender   Not marked Marked
##   Male           -5      5
##   Female          5     -5

corrplot(test$stdres, is.cor = FALSE)

Positive residuals (blue cells) show a positive association between the gender and being exposed to breathing smoke in a job, while negative residuals (red cells) show a negative association. It means that males are positively associated with being exposed to breathing smoke in a job and females are positively associated with not being exposed to breathing smoke in a job.

assocplot(t(cont_table), main = "Residuals and number of observations")

The cells for females which are not exposed to breathing smoke in a job and for males who are (black color) indicate that the observations are higher than it was expected. Other two cells (red color) show that the observed values are lower than expected ones. It means that there are much more men than women who are exposed to breathing smoke in a job and much more women that men who are not exposed to breathing smoke in a job.

Conclusion

We decline the H0-hypothesis and conclude that gender and being exposed to breathing smoke in a job are dependent variables. Men are more likely to be exposed to breathing smoke in a job, while women are not.

Hypotheses

hh0 <- as.data.frame(cbind("H0:", "The average number of cigarettes smoked per day by females is not different from that smoked by males."))
hh1 <- as.data.frame(cbind("H1:", "The average number of cigarettes smoked per day by females is different from that smoked by males."))
HHs <- as.data.frame(rbind(hh0, hh1))
colnames(HHs) <- c(" ", " ")

HHs %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F) %>%
  row_spec(1:2, background = "White")

Data Examination

We build a boxpolot.

ggplot(data2, aes(y = cgtsday, x = gndr)) + 
  geom_boxplot(fill = c("skyblue", "salmon")) +
  stat_summary(fun.y = mean, geom = "point", shape = 4, size = 4) +
  theme_classic() +
  labs(title = "Difference in amount of cigerettes smoked a day between males and females", x = "Gender", y = "Cigarettes smoked per day")

From the boxplot it is seen that in Slovenia men smoke more cigarettes per day than women.

We check the distribution of our data with numbers.

describeBy(data2, data2$gndr)

## 
##  Descriptive statistics by group 
## group: Male
##          vars   n mean   sd median trimmed mad min max range skew kurtosis
## gndr*       1 139  1.0 0.00      1     1.0 0.0   1   1     0  NaN      NaN
## jbexebs*    2 139  1.3 0.47      1     1.3 0.0   1   2     1 0.68    -1.55
## cgtsday     3 139 15.3 7.54     17    15.3 4.4   1  40    39 0.01    -0.19
##            se
## gndr*    0.00
## jbexebs* 0.04
## cgtsday  0.64
## -------------------------------------------------------- 
## group: Female
##          vars   n mean   sd median trimmed mad min max range skew kurtosis
## gndr*       1 148  2.0 0.00      2       2 0.0   2   2     0  NaN      NaN
## jbexebs*    2 148  1.1 0.29      1       1 0.0   1   2     1 2.74     5.56
## cgtsday     3 148 11.1 6.45     10      11 7.4   1  30    29 0.45    -0.53
##            se
## gndr*    0.00
## jbexebs* 0.02
## cgtsday  0.53

For male skew value is 0.01, so as should be for normal distribution. For female it is 0.45, so less normally but still less than 1.0. The negative kurtosis values show that the data has lighter tails than the standard distribution.

We build a histogram.

ggplot(data2, aes(x = cgtsday, fill = gndr), na.rm = TRUE) +
  geom_histogram(binwidth = 5, alpha = 0.75) +
  geom_density(alpha = 0.5) +
  labs(title = "Smoking habits in Slovenia",x = "Cigarettes smoked per day", y = "Density") +
  geom_vline(aes(xintercept = mean(data2$cgtsday, na.rm = TRUE), color='mean'), show.legend = TRUE, size = 1) +
  geom_vline(aes(xintercept = median(data2$cgtsday, na.rm = TRUE), color='median'), show.legend = TRUE, size = 1) +
  scale_fill_manual(values = c("skyblue", "salmon"), guide = FALSE) +
  scale_color_manual(values = c("hotpink4", "orangered")) +
  theme(legend.title = element_blank()) +
  facet_grid(data2$gndr)

Histogram shows that distribution in both groups is close to normal but still is not.

We filter our data by gender and build Q-Q plot.

m <- data2 %>%
  filter(gndr == 'Male')
f <- data2 %>%
  filter(gndr == 'Female')
par(mfrow = c(1,2))
qqnorm(m$cgtsday); qqline(m$cgtsday, col= "red", lty = 5, lwd = 2)
qqnorm(f$cgtsday); qqline(f$cgtsday, col = "red", lty = 5, lwd = 2)

Q-Q plots show that the distributions are not really normal.

T-Test

We apply T-test formula.

t.test(f$cgtsday, m$cgtsday, paired = F, var.equal = F)

## 
##  Welch Two Sample t-test
## 
## data:  f$cgtsday and m$cgtsday
## t = -5, df = 300, p-value = 8e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.8 -2.6
## sample estimates:
## mean of x mean of y 
##        11        15

On average men smoke 15 cigarettes a day whereas women smoke only 11. The t-statistic t(272.2) = - 5 (p-value < 0.001), so we can claim that there is a statistically significant difference between men and women in amount of cigarettes they smoke in a day.

Wilcox Test

We apply non-parametric test.

wilcox.test(x = f$cgtsday, y = m$cgtsday, paired = F)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  f$cgtsday and m$cgtsday
## W = 7000, p-value = 2e-06
## alternative hypothesis: true location shift is not equal to 0

The Wilcoxon rank sum W = 6985.5 (p-value < 0.001), which means that different genders smoke different amount of cigarettes in a day and this difference is statistically significant.

Conclusion

We reject the H0-hypothesis and conclude that the average number of cigarettes smoked per day by males is bigger than that smoked by females.

Marital Status

par(mar = c(5, 25, 2, 2))
barplot(table(data3$maritalb)/nrow(data3)*100, las = 2, horiz = T, col = "cornflowerblue")  

Given our analysis is aiming to understand whether marital status has an effect on weight we have decided to merge variable “Legally married”" with variable “In a legally registered civil union” and variable “Legally divorced/civil union dissolved” with variable “Legally separated”. For the purpose of this analysis they are the same - a civil union is a non religious marriage certificate (often between homosexual couples whereby the church, or the law, prohibits marriages). Slovenia has recognized partnerships (Slovene: “partnerska zveza”“) since 24 February 2017. These provide same-sex partners with all the legal rights of marriages, with the exception of joint adoption and in vitro fertilisation. Previously, Slovenia had recognized the more limited”registrirana partnerska skupnost" for same-sex couples since 23 July 2006, which gave same-sex partners access to one another’s pensions and property.

data3$maritalb1 <- ifelse(data3$maritalb == "Legally married", "Married",
                             ifelse(data3$maritalb == "In a legally registered civil union", "Married",
                                    ifelse(data3$maritalb == "Legally divorced/civil union dissolved", "Divorced",
                                           ifelse(data3$maritalb == "Legally separated", "Divorced",
                                                  ifelse(data3$maritalb == "Widowed/civil partner died", "Widowed", "Never married")
                                           )
                                    )
                             )
)

par(mar = c(3,10,0,3))
barplot(table(data3$maritalb1)/nrow(data3)*100, las = 2, horiz = T, col = "cornflowerblue")

Now we have four groups of comparable size, which is needed for ANOVA.

Weight by Marital Status

describeBy(data3$weight1, data3$maritalb1, mat = TRUE) %>% 
  select(Marital_Status = group1, N=n, Mean=mean, SD=sd, Median=median, Min=min, Max=max, 
                Skew=skew, Kurtosis=kurtosis, st.error = se) %>% 
  kable(align=c("lrrrrrrrr"), digits=2, row.names = FALSE,
        caption="Weight by Marital Status") %>% 
  kable_styling(bootstrap_options=c("bordered", "responsive","striped"), full_width = FALSE) 

The weight is normally distrubited within different marital statuses.

Then we create a boxplot.

ggplot(data3, aes(x = maritalb1, y = weight1)) +
  geom_boxplot(fill= "cornflowerblue") +
  theme_classic() +
  labs(x = "Marital Status", y = "Weight", title = "The Distrubution of Weight according to Marital Status") +
  theme(legend.position="none")

From the boxplot, we see that the weight is distributed rather normally across the marital statuses and that it is slightly lower among the never married respondents.

Homogenity of Variances

leveneTest(data3$weight1 ~ data3$maritalb1)  

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value Pr(>F)
## group    3    1.64   0.18
##       1180

P-value is 0.1783 therefore variances are equal. Thus, in ANOVA we should indicate var.equal = T.

F-test

aov(data3$weight1 ~ data3$maritalb1)

## Call:
##    aov(formula = data3$weight1 ~ data3$maritalb1)
## 
## Terms:
##                 data3$maritalb1 Residuals
## Sum of Squares             9220    277629
## Deg. of Freedom               3      1180
## 
## Residual standard error: 15
## Estimated effects may be unbalanced

aov.out <- aov(data3$weight1 ~ data3$maritalb1) 
summary(aov.out) 

##                   Df Sum Sq Mean Sq F value  Pr(>F)    
## data3$maritalb1    3   9220    3073    13.1 2.1e-08 ***
## Residuals       1180 277629     235                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F(3, 1180) = 13.062 and p-value of <0.001 indicate that the difference in weight across marital status is statistically significant.

Checking Residuals

Let’s check normality of our residuals.

plot(aov.out, 2) 

Residual are normal.

Post hoc Test

Next, we want to find out which of our variables create the significance. As our variances are equal we can run a Tukey test to find out.

par(mar = c(5, 10, 3, 1))
Tukey <- TukeyHSD(aov.out)
plot(Tukey, las = 2)  

As we can see from the graph, the pairwise tests that have significantly different means at a 95% confidence level are ‘Married’ and ‘Never married’ and ‘Widowed’ and ‘Married’. This fits the sociological theory that marrying and being widowed cause a change in weight.

Non-parametric Test

Non-parametric test for ANOVA is the Kruskal-Wallis test.
H0: Mean ranks of people from different marital satuses are the same.

data3$maritalb2 <- as.factor(data3$maritalb1)
kruskal.test(data3$weight1 ~ data3$maritalb2)  

## 
##  Kruskal-Wallis rank sum test
## 
## data:  data3$weight1 by data3$maritalb2
## Kruskal-Wallis chi-squared = 40, df = 3, p-value = 2e-08

With KW chi-square(3) = 38.916 and p-value < 0.001, we can state that the mean ranks of different marital statuses are not the same. This results confirm our results in the ANOVA test.

As Kruskal-Wallis test is significant, we run non-parametric post hoc test - Dunn?s test.

DunnTest(weight1 ~ maritalb2, data = data3) 

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                        mean.rank.diff    pval    
## Married-Divorced                  2.2  0.9634    
## Never married-Divorced         -130.3  0.0331 *  
## Widowed-Divorced                -84.5  0.3902    
## Never married-Married          -132.5 8.7e-09 ***
## Widowed-Married                 -86.7  0.0635 .  
## Widowed-Never married            45.8  0.4388    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Applying Dunn’s non-paramtric test, we can say that there is statistically significant difference between ‘Never Married - Divorced’, ‘Never married - Married’, and ‘Widowed-Married’, which actually proves the results of Kruskal-Wallis test.

Data Examination

data4 <- select(ESS, alcwkdy, alcwknd, gndr)
data4 <- na.omit(data4)

Firstly, we check if our variables are of the correct type.

str(data4)

## 'data.frame':    938 obs. of  3 variables:
##  $ alcwkdy: Factor w/ 80 levels "0","4","7","9",..: 15 1 2 1 29 9 32 19 11 19 ...
##  $ alcwknd: Factor w/ 87 levels "0","4","7","9",..: 12 5 2 11 38 9 5 32 18 25 ...
##  $ gndr   : Factor w/ 2 levels "Male","Female": 1 1 2 1 1 2 1 1 1 1 ...
##  - attr(*, "variable.labels")= Named chr  "Title of dataset" "ESS round" "Edition" "Production date" ...
##   ..- attr(*, "names")= chr  "name" "essround" "edition" "proddate" ...
##  - attr(*, "codepage")= int 65001
##  - attr(*, "na.action")= 'omit' Named int  1 3 10 11 12 13 15 16 17 18 ...
##   ..- attr(*, "names")= chr  "1" "3" "10" "11" ...

All the variables except for ‘gender’ should be numeric. Thus, we changed the type of these variables.

for (i in 1:2)
{
  data4[,i] <- as.numeric(as.character(data4[,i]))
}  

summary(data4)

##     alcwkdy       alcwknd        gndr    
##  Min.   :  0   Min.   :  0   Male  :466  
##  1st Qu.: 10   1st Qu.: 10   Female:472  
##  Median : 13   Median : 19               
##  Mean   : 22   Mean   : 28               
##  3rd Qu.: 26   3rd Qu.: 38               
##  Max.   :354   Max.   :286

Now, it looks fine!

We decided to present our description of chosen variables in form of the table.

Gender

ggplot(data4, aes(x = gndr)) + 
  geom_bar(aes(y = (..count..)/sum(..count..)), legend = TRUE, fill = c("skyblue", "salmon")) + 
  scale_y_continuous(labels=scales::percent) +
  ggtitle("Gender of the Respondent") + 
  xlab("Gender") + 
  ylab(" ")

As we can see, amount of males and females in our data is almost the same, so we do not need to apply set.seed() function in order to make it equal.

Alcohol Consumption

Let’s firstly build graphs of alcohol consumption on weekday and weekends.

ggplot(data4, aes(x = alcwkdy)) + 
  geom_density(fill = "gold1", colour = "goldenrod2", alpha = 0.6) + 
  ggtitle("Alcohol Consumption on a Weekday") +
  xlab("grams of alcohol") +
  scale_x_continuous(breaks = seq(0, 350, 25), limits=c(0, 354))

ggplot(data4, aes(x = alcwknd)) + 
  geom_density(fill = "#4271AE", colour = "#1F3552", alpha = 0.6) + 
  ggtitle("Alcohol Consumption on a Weekend") +
  xlab("grams of alcohol") +
  scale_x_continuous(breaks = seq(0, 250, 25), limits=c(0, 286))

Both distributions are not destributed normaly and both right-tailed. Let’s also build qqplot separated for both variables.

QQ Plot and T-Test

qqnorm(data4$alcwknd); qqline(data4$alcwknd, col= "red", lty = 5, lwd = 2)

qqnorm(data4$alcwkdy); qqline(data4$alcwkdy, col = "red", lty = 5, lwd = 2)

Distributions are not normal as the line is not straight.

t.test(data4$alcwkdy, data4$alcwknd, paired = F, var.equal = F)

## 
##  Welch Two Sample t-test
## 
## data:  data4$alcwkdy and data4$alcwknd
## t = -5, df = 2000, p-value = 5e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -8.8 -3.5
## sample estimates:
## mean of x mean of y 
##        22        28

Difference in the amount of alcohol last consumed on a weekday and weekend does has statistical significance, with a p value >0.001.

Consumption by Gender

ggplot(data4, aes(y = alcwkdy, x = gndr)) + 
  geom_boxplot(fill = c("skyblue", "salmon")) +
  stat_summary(fun.y = mean, geom = "point", shape = 4, size = 4) +
  theme_classic() +
  labs(title = "Weekday Drinking Between Males and Females", x = "Gender", y = "Amount of alcohol last drunk on a weekday")

ggplot(data4, aes(y = alcwknd, x = gndr)) + 
  geom_boxplot(fill = c("skyblue", "salmon")) +
  stat_summary(fun.y = mean, geom = "point", shape = 4, size = 4) +
  theme_classic() +
  labs(title = "Weekend Drinking Between Males and Females", x = "Gender", y = "Amount of alcohol last drunk on a weekend")

As it is seen from the boxplots males drink more than females of both weekends and weekdays.

Model 1

The F-statistic is significant (p-value: < 2.2e-16), therefore we can say that our model is better than no model at all. According to adjusted R-squared, which shows the explanatory power of a model and estimates 0.09, we can claim that the given model explains aproximately 9% of the variance. Here the intercept is the predicted value of grams of alcohol consumed at weekends, when amount of grams of alcohol consumed at weekdays is equal to 0, here it estimates 20.6 grams. Regression equation for the first model:

Predicted Grams Of Alcohol On Weekend Day = 20,6 + 0,33 * Grams Of Alcohol On Weekday

The interprtation:
With each increase in grams of alcohol drank on weekday by one, grams of alcohol drank on weekend day increases by 0.33.

plot(model1, which=1)  

The graph ‘Residuals vs Fitted’ shows that residuals are distributed not quite linearly because the red line is not straight enough and that the residuals are also not homoscedastic as the points are not equally distributed around the 0 line.

Model 2

F-statistic is statistically significant as the p-value < 2.2e-16 and to be significant it needs to be less than 0.05. Thus, we can conclude that our model is better than the one that predicts the amount of alcohol that people drink on weekends with mean value. From the adjusted R-squared we can conclude that our second model explains 13% of variation in grams of alcohol people drink on weekends. The intercept is the predicted value of grams of alcohol drunk on weekends when the amount of grams of alcohol drank on weekdays is equal to 0 and a person is a male, and it is equal to 27,89 grams. So, our regression equation looks like this:

Predicted Grams Of Alcohol On Weekend Day = 27,89 + 0,29 * Grams Of Alcohol On Weekday - 12,53 * Female

The interpretation:
With each increase in grams of alcohol drank on weekday by one, grams of alcohol drank on weekend day increases by 0.29 and if a respondent is a female, the amount of grams of alcohol drank on weekend day decreases by 12.53.

plot(model2, which=1)

The graph illustrates that residuals are distributed almost linearly because the red line is very close to the straight one and that the residuals are not homoscedastic as the points are not equally distributed around the 0 line.

Model 3

F-statistic is statistically significant as the p-value < 2.2e-16 and to be significant it needs to be less than 0.05. From the adjusted R-squared we can conclude that the third model explains almost 15% of variation in grams of alcohol people drink on weekends.

Predicted Grams Of Alcohol On Weekend Day (Female) = 24,64 + 0,41 * Grams Of Alcohol On Weekday - 6,98 * Female - 0.26 * Grams Of Alcohol On Weekday * Female

Predicted Grams Of Alcohol On Weekend Day (Male) = 24,64 + 0,41 * Grams Of Alcohol On Weekday

To explain our interaction effect: alcohol consumption on weekends actually depends not solely on weekdays consumption OR just on your gender, but on the interaction of these two variables. So the association is positive and we see the more you drink on the weekdays the more you drink on the weekend, but the rate of increase is larger for males than females. However, confusingly and unexpectedly, the amount consumed is higher during the week than during the weekend.

To summarise this is plain english: for females, with each additional gram of alcohol on weekday, weekend consumption of alcohol increases on 0.41 gram in comparison with the baseline group (males who don’t drink on weekdays*, in this instance).

For example:

A female, who consumes 300 grams of alcohol during week days is predicted to consume on weekend days the following amount of alcohol (in grams):
\(24.64 + 0.41 * 300 - 6.98 - 0.26 * 300 = 62.66\)
62.66 predicted grams of alcohol consumed by a female on weekend days if she consumed 300 grams of alcohol on week days.

A male, who consumes 300 grams of alcohol during week days is predicted to consume on weekend daysthe following amount of alcohol (in grams):
\(24.64 + 0.41 * 300 = 147.64\)
147.64 predicted grams of alcohol consumed by a male on weekend days if he consumed 300 grams of alcohol on week days.

plot(model3, which = 1)  

Here our plot suggests that our residuals (measured by the red line) are most closely fitted to 0.

Calculating Accuracy

Having made our models and examined the details of our various models we should try and deduce which model is the best fit. To do this we will use both the ANOVA function, to compare the statistical difference in the means, and then look at quantifying how well the model fits the data. Firstly, lets run ANOVA.

anova(model1, model2, model3)

## Analysis of Variance Table
## 
## Model 1: alcwknd ~ alcwkdy
## Model 2: alcwknd ~ alcwkdy + gndr
## Model 3: alcwknd ~ alcwkdy * gndr
##   Res.Df    RSS Df Sum of Sq    F  Pr(>F)    
## 1    936 777588                              
## 2    935 742279  1     35309 45.2 3.2e-11 ***
## 3    934 730248  1     12031 15.4 9.4e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We have a p value of <0.001. This suggests that there is a significant difference between the mean values of our ANOVAs. In other words - our third model, which has an interaction effect between gender and alcohol consumption is better than our first. But by how much?

One way of understanding what this actually means - and by how much our models differ - is to run a comparison on the sum of square residuals and calculate the square root of the mean of residuals^2. The first gives us an overall figure to compare how much our model is out by and the second gives the average difference between our model and the data point.

rgrs <- list(model1, model2, model3)

ssq_function <- function(lm) {
  get_regression_points(lm) %>%
    mutate(sq_residuals = residual^2) %>%
    summarize(sum_sq_residuals = sum(sq_residuals), rmse = sqrt(mean(sq_residuals)))
}

for(i in 1:3) {
  print(ssq_function(rgrs[[i]]))
}

## # A tibble: 1 x 2
##   sum_sq_residuals  rmse
##              <dbl> <dbl>
## 1          777586.  28.8
## # A tibble: 1 x 2
##   sum_sq_residuals  rmse
##              <dbl> <dbl>
## 1          742279.  28.1
## # A tibble: 1 x 2
##   sum_sq_residuals  rmse
##              <dbl> <dbl>
## 1          730248.  27.9

As we see, the sum of square residuals is lowest for our third model.

Train-Test Split

Alternatively, we can cross validate our results by running a train/test split.

To test the accuracy of our model we’re going to split up the data into two parts (based on a roughly 80/20 split, because of the relatively small size of the data set we have chosen to improve the accuract of our model by increasing the size of the training set). We then create a linear model based on the training data. We then compare the predicted results based on our training results to the remaining 20% of the dataset, calculating the square root of the average residual to get an idea of how accurate our model is.

set.seed(76)

data_shuffled <- data4 %>% 
sample_frac(size = 1, replace = FALSE)

# Train/test split
train <- data_shuffled %>%
slice(1:800)
test <- data_shuffled %>%
slice(801:983)

# Fit model to training set
train_model_2 <- lm(alcwknd ~ alcwkdy * gndr, data=train)

We’ve created the linear model based on 80% of the data. Now we’re going to compare the residuals between the remaining 20% of the data and the values our model has predicted (in accordance with the train data).

get_regression_points(train_model_2, newdata = test)

## # A tibble: 138 x 6
##       ID alcwknd alcwkdy gndr   alcwknd_hat residual
##    <int>   <dbl>   <dbl> <fct>        <dbl>    <dbl>
##  1     1      10       0 Male          25.1   -15.1 
##  2     2      48      10 Female        19.1    28.9 
##  3     3      26      29 Female        21.5     4.55
##  4     4      60      60 Male          48.5    11.5 
##  5     5      14       4 Female        18.4    -4.39
##  6     6      19      19 Female        20.2    -1.23
##  7     7      46      33 Male          38       8   
##  8     8      23      20 Male          32.9    -9.92
##  9     9      10      10 Female        19.1    -9.13
## 10    10      29      29 Male          36.4    -7.44
## # ... with 128 more rows

get_regression_points(train_model_2, newdata = test) %>% 
mutate(sq_residuals = residual^2) %>%
summarize(rmse = sqrt(mean(sq_residuals)))

## # A tibble: 1 x 1
##    rmse
##   <dbl>
## 1  31.3

Making Real Life Predictions

What, though, does all of this mean?

In accordance with Bernardi et al’s best practice for reporting statistical significance we’re going to demonstrate the difference in our models.

An (imaginary) man comes to us and tells us he’s worried about how much his father drinks. We enquire how much his father usually drinks on a weekday. The man tells us he drinks 25g of alcohol. We use our model to predict how much he is likely to be drinking on a weekend day (the man, owing to working on the weekend, is unable to closely monitor his father, but the empty vodka bottles he finds the following morning are worrying…) and thus whether this man has any cause to be concerned.

model1_pred <- data.frame("gndr" = 'Male', "alcwkdy" = 25)
augment(model1, newdata = model1_pred, type.predict = "response")

## # A tibble: 1 x 4
##   gndr  alcwkdy .fitted .se.fit
##   <fct>   <dbl>   <dbl>   <dbl>
## 1 Male       25    28.9   0.948

196g of alcohol is the recommended weekly allowance according to US health guidelines. So far we know he drinks 100g of alcohol, on average, during the week. Based on our first (more inaccurate) model we would conclude that he drinks a further 87g over the weekend (3 * 29). This totals at 187g a week - under the limit.

model3_pred <- data.frame("gndr" = 'Male', "alcwkdy" = 25)
augment(model3, newdata = model3_pred, type.predict = "response")

## # A tibble: 1 x 4
##   gndr  alcwkdy .fitted .se.fit
##   <fct>   <dbl>   <dbl>   <dbl>
## 1 Male       25    34.8    1.30

However, when we run the same data through our more accurate third model we reach a different conclusion. Instead, we predict he drinks 105g (3 * 35) of alcohol over the weekend. This means his father drinks 205g, over the recommended limit.

In this instance, our more accurate model has allowed us to determine whether someone’s lifestyle has become harmful or not.

Finally, lets compare how accurate our prediction for this man has been compared to what our real data tells us. Lets wrangle our data to identify men who drink 25g of alcohol and what their average alcohol consumption is.

man <- data %>%
  filter(gndr == 'Male') %>%
  filter(alcwkdy >= 24) %>%
  filter(alcwkdy <= 26) %>%
  summarize("Number of observations" = n(), "Grams of alcohol drunk on the weekend" = mean(alcwknd))
man

##   Number of observations Grams of alcohol drunk on the weekend
## 1                     13                                    36

Comparing the real data to our prediction shows that the amount of alcohol drunk on the weekend was actually 35.6, compared to our models prediction of 34.8. Pretty accurate! Whilst our small sample size is only 13 it nonetheless shows the accuracy of using linear models for predictions.