report_draft
2024-03-29
SOC222 Research Report Data for Marina Wilson (1009109901), Aliyah Fatima Walli (1008802770), and MeryemKarayunusoglu. In this report, we will be examining how education (Ba vs no BA) and religiouss belonging relate to happiness in Turkey.
Division of Labour Statement: All coding done by Marina Wilson (1009109901). Division of written content can be found in the appendix section of the report.
## Installing necessary packages.
packages <- c("tidyverse", "infer", "fst", "broom", "modelsummary", "effects", "survey", "MASS", "aod", "interactions", "kableExtra", "flextable", "scales") # add any you need here
# Install packages if they aren't installed already
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)
# Load the packages
lapply(packages, library, character.only = TRUE)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors## Warning: package 'infer' was built under R version 4.3.3## Warning: package 'effects' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3## lattice theme set by effectsTheme()
## See ?effectsTheme for details.## Warning: package 'survey' was built under R version 4.3.3## Loading required package: grid
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Loading required package: survival
##
## Attaching package: 'survey'
##
## The following object is masked from 'package:graphics':
##
## dotchart
##
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select## Warning: package 'aod' was built under R version 4.3.3##
## Attaching package: 'aod'
##
## The following object is masked from 'package:survival':
##
## rats## Warning: package 'interactions' was built under R version 4.3.3##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows
##
##
## Attaching package: 'flextable'
##
## The following objects are masked from 'package:kableExtra':
##
## as_image, footnote
##
## The following object is masked from 'package:purrr':
##
## compose
##
##
## Attaching package: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factor## [[1]]
## [1] "lubridate" "forcats" "stringr" "dplyr" "purrr" "readr"
## [7] "tidyr" "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [13] "grDevices" "utils" "datasets" "methods" "base"
##
## [[2]]
## [1] "infer" "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [7] "readr" "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [13] "graphics" "grDevices" "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "fst" "infer" "lubridate" "forcats" "stringr" "dplyr"
## [7] "purrr" "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [13] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [19] "base"
##
## [[4]]
## [1] "broom" "fst" "infer" "lubridate" "forcats" "stringr"
## [7] "dplyr" "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [13] "tidyverse" "stats" "graphics" "grDevices" "utils" "datasets"
## [19] "methods" "base"
##
## [[5]]
## [1] "modelsummary" "broom" "fst" "infer" "lubridate"
## [6] "forcats" "stringr" "dplyr" "purrr" "readr"
## [11] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [16] "graphics" "grDevices" "utils" "datasets" "methods"
## [21] "base"
##
## [[6]]
## [1] "effects" "carData" "modelsummary" "broom" "fst"
## [6] "infer" "lubridate" "forcats" "stringr" "dplyr"
## [11] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [16] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [21] "datasets" "methods" "base"
##
## [[7]]
## [1] "survey" "survival" "Matrix" "grid" "effects"
## [6] "carData" "modelsummary" "broom" "fst" "infer"
## [11] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [16] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base"
##
## [[8]]
## [1] "MASS" "survey" "survival" "Matrix" "grid"
## [6] "effects" "carData" "modelsummary" "broom" "fst"
## [11] "infer" "lubridate" "forcats" "stringr" "dplyr"
## [16] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [21] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [26] "datasets" "methods" "base"
##
## [[9]]
## [1] "aod" "MASS" "survey" "survival" "Matrix"
## [6] "grid" "effects" "carData" "modelsummary" "broom"
## [11] "fst" "infer" "lubridate" "forcats" "stringr"
## [16] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [21] "ggplot2" "tidyverse" "stats" "graphics" "grDevices"
## [26] "utils" "datasets" "methods" "base"
##
## [[10]]
## [1] "interactions" "aod" "MASS" "survey" "survival"
## [6] "Matrix" "grid" "effects" "carData" "modelsummary"
## [11] "broom" "fst" "infer" "lubridate" "forcats"
## [16] "stringr" "dplyr" "purrr" "readr" "tidyr"
## [21] "tibble" "ggplot2" "tidyverse" "stats" "graphics"
## [26] "grDevices" "utils" "datasets" "methods" "base"
##
## [[11]]
## [1] "kableExtra" "interactions" "aod" "MASS" "survey"
## [6] "survival" "Matrix" "grid" "effects" "carData"
## [11] "modelsummary" "broom" "fst" "infer" "lubridate"
## [16] "forcats" "stringr" "dplyr" "purrr" "readr"
## [21] "tidyr" "tibble" "ggplot2" "tidyverse" "stats"
## [26] "graphics" "grDevices" "utils" "datasets" "methods"
## [31] "base"
##
## [[12]]
## [1] "flextable" "kableExtra" "interactions" "aod" "MASS"
## [6] "survey" "survival" "Matrix" "grid" "effects"
## [11] "carData" "modelsummary" "broom" "fst" "infer"
## [16] "lubridate" "forcats" "stringr" "dplyr" "purrr"
## [21] "readr" "tidyr" "tibble" "ggplot2" "tidyverse"
## [26] "stats" "graphics" "grDevices" "utils" "datasets"
## [31] "methods" "base"
##
## [[13]]
## [1] "scales" "flextable" "kableExtra" "interactions" "aod"
## [6] "MASS" "survey" "survival" "Matrix" "grid"
## [11] "effects" "carData" "modelsummary" "broom" "fst"
## [16] "infer" "lubridate" "forcats" "stringr" "dplyr"
## [21] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [26] "tidyverse" "stats" "graphics" "grDevices" "utils"
## [31] "datasets" "methods" "base"## Reading in the Turkey data and creating a dataframe for it.
turkey_data <- read.fst("turkey_data.fst")
df <- turkey_dataSection 1: Determining Descriptive Statistics of the Dataset.
## First, I must examine and clean all variables I will be working with; starting with happiness.
turkey_happy <- df %>%
filter(cntry == "TR") %>%
dplyr:: select(happy)
turkey_happy$y <- turkey_happy$happy
table(turkey_happy$y)##
## 0 1 2 3 4 5 6 7 8 9 10 88
## 202 117 146 226 291 887 432 530 540 296 517 88## 88 must be removed from the dataset, or results will be skewed.
turkey_happy$y[turkey_happy$y %in% 88] <- NA
table(turkey_happy$y)##
## 0 1 2 3 4 5 6 7 8 9 10
## 202 117 146 226 291 887 432 530 540 296 517## Next, I will generate a table of descriptive statistics for the variable "how happy are you" using the datasummary_skim function.
tablehap <- datasummary_skim ((turkey_happy$y), title = "Table 1. Descriptive Statistics of the -How Happy Are You- Variable", output = "flextable")## Warning: The histogram argument is only supported for (a) output types "default",
## "html", "kableExtra", or "gt"; (b) writing to file paths with extensions
## ".html", ".jpg", or ".png"; and (c) Rmarkdown, knitr or Quarto documents
## compiled to PDF (via kableExtra) or HTML (via kableExtra or gt). Use
## `histogram=FALSE` to silence this warning.tablehap
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
data | 12 | 2 | 6.0 | 2.7 | 0.0 | 6.0 | 10.0 |
## Now, I will repeat the same steps for our second variable, "belonging to a particular religion or denomination".
turkey_rlgblg <- df %>%
filter(cntry == "TR") %>%
dplyr:: select(rlgblg)
turkey_religion <- turkey_rlgblg$rlgblg
table(turkey_religion)## turkey_religion
## 1 2 8
## 4135 127 10## Cleaning religion to remove NA's (8).
turkey_religion[turkey_religion %in% 8] <- NA
table(turkey_religion)## turkey_religion
## 1 2
## 4135 127## Now generating a table of descriptive statistics for the variable "belonging to a particular religion or denomination".
tablerlg <- datasummary_skim ((turkey_religion), title = "Table 2. Descriptive Statistics of the -Belonging to a Particular Religion or Denomination- Variable", output = "flextable")## Warning: The histogram argument is only supported for (a) output types "default",
## "html", "kableExtra", or "gt"; (b) writing to file paths with extensions
## ".html", ".jpg", or ".png"; and (c) Rmarkdown, knitr or Quarto documents
## compiled to PDF (via kableExtra) or HTML (via kableExtra or gt). Use
## `histogram=FALSE` to silence this warning.tablerlg
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
data | 3 | 0 | 1.0 | 0.2 | 1.0 | 1.0 | 2.0 |
## Finally, repeating the same steps for education level (BA or no BA). Starting with checking the data. Using edulvla rather than edulvlb, because no results come up for Turkey with edulvlb.
table(df$edulvla)##
## 1 2 3 5 55 88 99
## 2481 680 787 314 8 1 1## Now cleaning the variables edulvlb and edulvla.
df <- df %>%
# Set specific values in 'edulvla' and 'edulvlb' to NA
mutate(
# Replace 55, 77, and 88 in 'edulvla' with NA
edulvla = na_if(edulvla, 55) %>%
na_if(., 88) %>%
na_if(., 99),
# Replace 5555, 7777, and 8888 in 'edulvlb' with NA
edulvlb = na_if(edulvlb, 5555) %>%
na_if(., 7777) %>%
na_if(., 8888)
)
## Checking the data:
table(df$edulvla)##
## 1 2 3 5
## 2481 680 787 314## Now generating a table of descriptive statistics for the variable "highest level of education".
tableedu <- datasummary_skim ((df$edulvla), title = "Table 3. Descriptive Statistics of the -Highest Level of Education- Variable", output = "flextable")## Warning: The histogram argument is only supported for (a) output types "default",
## "html", "kableExtra", or "gt"; (b) writing to file paths with extensions
## ".html", ".jpg", or ".png"; and (c) Rmarkdown, knitr or Quarto documents
## compiled to PDF (via kableExtra) or HTML (via kableExtra or gt). Use
## `histogram=FALSE` to silence this warning.tableedu
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
data | 5 | 0 | 1.8 | 1.2 | 1.0 | 1.0 | 5.0 |
## Now creating a table of descriptive statistics for all three variables (both predictors and outcome) for comparison reasons. Will be used as Table 1 in the report.
turkey_clean <- df %>%
mutate(
happy = ifelse(happy %in% c(88), NA, happy),
rlgblg = ifelse(rlgblg %in% c(8), NA, rlgblg),
edulvla = ifelse(edulvla %in% c(55, 77, 88), NA, edulvla),
)
turkey_clean <- turkey_clean %>% filter(!is.na(happy), !is.na(rlgblg), !is.na(edulvla))
table1 <- datasummary_skim((turkey_clean %>% dplyr::select(happy, rlgblg, edulvla)),title = "Table 1. Descriptive Statistics for the Variables of Happiness, Religious Belonging, and Education Level in Turkey", output = "flextable")## Warning: The histogram argument is only supported for (a) output types "default",
## "html", "kableExtra", or "gt"; (b) writing to file paths with extensions
## ".html", ".jpg", or ".png"; and (c) Rmarkdown, knitr or Quarto documents
## compiled to PDF (via kableExtra) or HTML (via kableExtra or gt). Use
## `histogram=FALSE` to silence this warning.table1
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
happy | 11 | 0 | 6.0 | 2.7 | 0.0 | 6.0 | 10.0 |
rlgblg | 2 | 0 | 1.0 | 0.2 | 1.0 | 1.0 | 2.0 |
edulvla | 4 | 0 | 1.8 | 1.2 | 1.0 | 1.0 | 5.0 |
## Renaming variables so that the table can be more easily interpreted.
turkey_clean <- turkey_clean %>%
mutate(
Happiness = happy,
Religion = rlgblg,
Education = edulvla
)
table1a <- datasummary_skim((turkey_clean %>% dplyr::select(Happiness, Religion, Education)),title = "Table 1. Descriptive Statistics for the Variables of Happiness, Religious Belonging, and Education Level in Turkey", output = "flextable")## Warning: The histogram argument is only supported for (a) output types "default",
## "html", "kableExtra", or "gt"; (b) writing to file paths with extensions
## ".html", ".jpg", or ".png"; and (c) Rmarkdown, knitr or Quarto documents
## compiled to PDF (via kableExtra) or HTML (via kableExtra or gt). Use
## `histogram=FALSE` to silence this warning.table1a
| Unique (#) | Missing (%) | Mean | SD | Min | Median | Max |
|---|---|---|---|---|---|---|---|
Happiness | 11 | 0 | 6.0 | 2.7 | 0.0 | 6.0 | 10.0 |
Religion | 2 | 0 | 1.0 | 0.2 | 1.0 | 1.0 | 2.0 |
Education | 4 | 0 | 1.8 | 1.2 | 1.0 | 1.0 | 5.0 |
## Saving the table as a pdf file. Used code below, but added # afterwards to avoid issues/halting when knitting the document to Rpubs.
##flextable::save_as_docx(table1a, path = "Table1New.docx",
##width = 5.0, height = 7.0)
##set_flextable_defaults(fonts_ignore=TRUE)
##print(table1, preview = "pdf")Section 2: Using the infer package to generate the null distribution (with confidence interval) and the calculated test statistic.
## First, I have to recode our three variables.
## Recoding "Belonging to a Particular Religion or Denomination" (rlgblg) as "Yes" or "No".
df <- df %>%
mutate(religion = case_when(
rlgblg == 2 ~ "No",
rlgblg == 1 ~ "Yes",
rlgblg %in% c(8) ~ NA_character_,
TRUE ~ as.character(rlgblg)
))
# check
table(df$religion)##
## No Yes
## 127 4135## Recoding "Highest Level of Education" (edulvla) to "BA" or "no BA".
df <- df %>%
mutate (
educ.ba = case_when(
essround < 5 & edulvla == 5 ~ "BA or more",
essround >= 5 & edulvlb > 600 ~ "BA or more",
TRUE ~ "No BA"
)
) %>%
# Ensure all NAs are handled uniformly
mutate(
edulvla = ifelse(edulvla %in% c(55, 77, 88, 99), NA_integer_, edulvla),
edulvlb = ifelse(edulvlb %in% c(5555, 7777, 8888), NA_integer_, edulvlb),
educ.ba = as.character(educ.ba) # Force educ.ba to be treated as character
)
table(df$educ.ba)##
## BA or more No BA
## 314 3958## Ensuring happy is cleaned. Will not be recoding it, as the numerical scale from 0-10 works well to describe happiness level.
df <- df %>%
mutate(
happiness = ifelse(happy %in% c(88, 77), NA, happy)
)
table(df$happiness)##
## 0 1 2 3 4 5 6 7 8 9 10
## 202 117 146 226 291 887 432 530 540 296 517##Step 1: Calculating the test statistic of the sample (using education as our explanatory and happiness as our response, since education level is our "main" predictor). I decided to do an two-sample T test, since the explanatory variable (education) is made up of two categorical variables, and the response (happiness) is made up of one numeric variable. Based on the reference chart in lecture 6, a two-sided T test is the best fit for these variables.
test_stat <- df %>%
specify(explanatory = educ.ba,
response = happiness) %>%
hypothesize(null = "independence") %>%
calculate(stat = "t") ## Warning: Removed 88 rows containing missing values.## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "BA or more" - "No BA", or divided in the order "BA or more" / "No BA"
## for ratio-based statistics. To specify this order yourself, supply `order =
## c("BA or more", "No BA")` to the calculate() function.print(test_stat$stat)## t
## 1.387274## Step 2: Simulating the null distribution.
null_distribution <- df %>%
specify(explanatory = educ.ba,
response = happiness) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "t")## Warning: Removed 88 rows containing missing values.## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "BA or more" - "No BA", or divided in the order "BA or more" / "No BA"
## for ratio-based statistics. To specify this order yourself, supply `order =
## c("BA or more", "No BA")` to the calculate() function.## Step 3: Calculating the p-value of the sample. Kept direction as "two-sided" since this aligns with the t-test on the reference chart (seen in lecture slideshow 6).
p_val <- null_distribution %>%
get_pvalue(obs_stat = test_stat, direction = "two-sided")
p_val## # A tibble: 1 × 1
## p_value
## <dbl>
## 1 0.17## Step 4: Visualization of the null distribution with confidence intervals.
conf_int <- null_distribution %>%
get_confidence_interval(level = 0.95, type = "percentile")
null_distribution %>%
visualize() +
shade_p_value(obs_stat = test_stat, direction = "two-sided") +
shade_confidence_interval(endpoints = conf_int)
null_distribution## Response: happiness (numeric)
## Explanatory: educ.ba (factor)
## Null Hypothesis: independence
## # A tibble: 1,000 × 2
## replicate stat
## <int> <dbl>
## 1 1 -0.441
## 2 2 1.34
## 3 3 -0.873
## 4 4 -0.657
## 5 5 0.0434
## 6 6 -1.50
## 7 7 0.488
## 8 8 -0.902
## 9 9 2.08
## 10 10 -0.240
## # ℹ 990 more rows## A p-value of 0.172 is extremely high, and does not fall below the significance level of p<0.05 (which would allow us to reject the null). Furthermore, the observed test statistic (the red line) is present within the null distribution, and the confidence interval falls within the null distribution as well. These results indicate that under frequentist hypothesis testing, we DO NOT reject the null hypothesis.Section 3: Regression Modelling with Single, Multivariate, and Interaction Models.
## First, the data must be weighted to ensure only the most accurate estimates are included.
df$weight <- df$dweight * df$pweight
survey_design <- svydesign(ids = ~1, data = df, weights = ~weight)## SLR Model: Generated a regression model using our single main predictor, education level (BA vs no BA).
model1 <- lm(happiness ~ educ.ba, data = df, weights = weight)
modelsummary(
list(model1),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "SLR Model: Education Level & Happiness",
statistic = NULL)| (1) | |
|---|---|
| (Intercept) | 6.2 (0.2)*** |
| educ.baNo BA | −0.2 (0.2) |
| Num.Obs. | 4184 |
| R2 | 0.000 |
| R2 Adj. | 0.000 |
| AIC | 20670.2 |
| BIC | 20689.2 |
| Log.Lik. | −10332.086 |
| RMSE | 2.69 |
## A coefficient of -0.2 for educ.ba (No BA) means that on average, Turks with no BA score 0.2 points lower on the "How Happy Are You" scale than Turks with a BA. However, there are no asterisks, which indicate there is a high probability of results being due to chance. ## MLR Model: Generated a regression model with both predictors (education level and religious belonging).
model2 <- lm(happiness ~ educ.ba + religion, data = df, weights = weight)
modelsummary(
list(model2),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "MLR Model: Education Level, Religious Belonging, & Happiness",
statistic = NULL)| (1) | |
|---|---|
| (Intercept) | 6.0 (0.3)*** |
| educ.baNo BA | −0.2 (0.2) |
| religionYes | 0.3 (0.3) |
| Num.Obs. | 4175 |
| R2 | 0.001 |
| R2 Adj. | 0.000 |
| AIC | 20630.9 |
| BIC | 20656.3 |
| Log.Lik. | −10311.470 |
| RMSE | 2.69 |
## religionYES indicates that Turks who report belonging to a particular region or denomination score 0.3 points higher on the happiness scale on average. Again, no asterisks, so there is a high probability this result may be due to chance.## MLR with Interaction Model: Generated a regression model with both predictors, plus interaction.
model3 <- lm(happiness ~ educ.ba + religion + educ.ba*religion, data = df, weights = weight)
modelsummary(
list(model3),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "MLR Interaction Model: Education Level, Religious Belonging, & Happiness",
statistic = NULL)| (1) | |
|---|---|
| (Intercept) | 5.1 (0.6)*** |
| educ.baNo BA | 0.9 (0.6) |
| religionYes | 1.3 (0.6)* |
| educ.baNo BA × religionYes | −1.2 (0.6)+ |
| Num.Obs. | 4175 |
| R2 | 0.002 |
| R2 Adj. | 0.001 |
| AIC | 20629.4 |
| BIC | 20661.1 |
| Log.Lik. | −10309.697 |
| RMSE | 2.69 |
## The score of -1.2 for NO BA x religionYES indicates that having no BA + having a religious affiliation, on average, results in one scoring 1.2 points lower on the happiness scale. Again, no asterisks. ## Final Step: Generating a table including all models for comparison and analysis purposes. Will be used in the research report.
modelsummary(
list(model1, model2, model3),
fmt = 1,
estimate = c( "{estimate} ({std.error}){stars}"),
title = "Linear Regression Models: SLR, MLR, & MLR Interaction for Education Level, Religious Belonging, & Happiness",
statistic = NULL)| (1) | (2) | (3) | |
|---|---|---|---|
| (Intercept) | 6.2 (0.2)*** | 6.0 (0.3)*** | 5.1 (0.6)*** |
| educ.baNo BA | −0.2 (0.2) | −0.2 (0.2) | 0.9 (0.6) |
| religionYes | 0.3 (0.3) | 1.3 (0.6)* | |
| educ.baNo BA × religionYes | −1.2 (0.6)+ | ||
| Num.Obs. | 4184 | 4175 | 4175 |
| R2 | 0.000 | 0.001 | 0.002 |
| R2 Adj. | 0.000 | 0.000 | 0.001 |
| AIC | 20670.2 | 20630.9 | 20629.4 |
| BIC | 20689.2 | 20656.3 | 20661.1 |
| Log.Lik. | −10332.086 | −10311.470 | −10309.697 |
| RMSE | 2.69 | 2.69 | 2.69 |
Section 4: Providing a Third Visual.
## For this section, I decided to generate a interaction plot to get a better look at how/if religion and educational level interact in Turkey.
interaction_plot <- effect("educ.ba*religion", model3, na.rm=TRUE)
plot(interaction_plot,
main="Interaction Effect",
xlab="Religiousness",
ylab="Happiness Scale")
interaction_plot##
## educ.ba*religion effect
## religion
## educ.ba No Yes
## BA or more 5.067609 6.350135
## No BA 5.950904 6.010905