Week 9 workbook and solutions

Joseph Bulbulia https://josephbulbulia.netlify.app (Victoria University of Wellington)https://www.wgtn.ac.nz
05-13-2021

Import the jittered NZAVS dataset

Show code 
### Libraries
library("tidyverse")
library("patchwork")
library("lubridate")
library("kableExtra")
library("gtsummary")
library("lubridate")
library("equatiomatic")
library("ggdag")
library("brms")
library("rstan")
library("rstanarm")
library("bayesplot")
library("easystats")
library("kableExtra")
library("broom")
library("tidybayes")
library("bmlm")
# if (!require(tidyLPA)) {
#   install.packages("tidyLPA")
# }
# rstan options
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores ())
theme_set(theme_classic())

# read data
# read data


nz_0 <- as.data.frame(readr::read_csv2(
  url(
    "https://raw.githubusercontent.com/go-bayes/psych-447/main/data/nzj.csv"
  )
))



# to relevel kessler 6 variables
f <-
  c(
    "None Of The Time",
    "A Little Of The Time",
    "Some Of The Time",
    "Most Of The Time",
    "All Of The Time"
  )



# get data into shape
nz_cr <- nz_0 %>%
  dplyr::mutate_if(is.character, factor) %>%
  select(
    -c(
      SWB.Kessler01,
      SWB.Kessler02,
      SWB.Kessler03,
      SWB.Kessler04,
      SWB.Kessler05,
      SWB.Kessler06
    )
  ) %>%
  dplyr::mutate(Wave = as.factor(Wave)) %>%
  dplyr::mutate(FeelHopeless = forcats::fct_relevel(FeelHopeless, f)) %>%
  dplyr::mutate(FeelDepressed = forcats::fct_relevel(FeelDepressed, f)) %>%
  dplyr::mutate(FeelRestless = forcats::fct_relevel(FeelRestless, f)) %>%
  dplyr::mutate(EverythingIsEffort = forcats::fct_relevel(EverythingIsEffort, f)) %>%
  dplyr::mutate(FeelWorthless = forcats::fct_relevel(FeelWorthless, f)) %>%
  dplyr::mutate(FeelNervous = forcats::fct_relevel(FeelNervous, f)) %>%
  dplyr::mutate(Wave = as.factor(Wave)) %>%
  dplyr::mutate(male_id = as.factor(Male)) %>%
  dplyr::mutate(date = make_date(year = 2009, month = 6, day = 30) + TSCORE) %>%
  dplyr::mutate(
    FeelWorthless_int = as.integer(FeelWorthless),
    FeelNervous_int =  as.integer(FeelNervous),
    FeelHopeless_int =  as.integer(FeelHopeless),
    EverythingIsEffort_int =  as.integer(EverythingIsEffort),
    FeelRestless_int =  as.integer(FeelRestless),
    FeelDepressed_int =  as.integer(FeelDepressed),
    HLTH.Fatigue_int = as.integer(HLTH.Fatigue + 1)
  ) %>%
  dplyr::mutate(yearS = TSCORE - min(TSCORE, na.rm = TRUE)) %>%
  dplyr::mutate(KESSLER6sum = as.integer(KESSLER6sum))

## 2018 only data
nz2018 <- nz_cr %>%
  dplyr::filter(Wave == 2018)

## 2019 only timeline
# Relevel covid timeline 2019
ord_dates_class_2019_only <- c("PreCOVID",
                               "JanFeb",
                               "EarlyMarch",
                               "Lockdown",
                               "PostLockdown")

nz2019 <- nz_cr %>%
  dplyr::filter(YearMeasured == 1) %>%
  dplyr::filter(Wave == 2019) %>%
  dplyr::group_by(Id) %>%
  dplyr::ungroup(Id) %>%
  dplyr::mutate(Covid_Timeline_cr =
                  as.factor(ifelse(
                    TSCORE %in% 3896:3921,
                    # feb 29 - march 25th
                    "EarlyMarch",
                    ifelse(
                      TSCORE %in% 3922:3954,
                      "Lockdown",
                      #march 26- Mon 27 April 2020
                      ifelse(
                        TSCORE > 3954,
                        # after april 27th 20202
                        "PostLockdown",
                        ifelse(TSCORE %in% 3842:3895,
                               # jan 6 to feb 28
                               "JanFeb",
                               "PreCOVID")
                      )
                    )
                  ))) %>%
  dplyr::mutate(Covid_Timeline_cr = forcats::fct_relevel(Covid_Timeline_cr, ord_dates_class_2019_only))



## 2018/2019 data set  
ord_dates_class <- c("Baseline",
                     "PreCOVID",
                     "JanFeb",
                     "EarlyMarch",
                     "Lockdown",
                     "PostLockdown")

nzl <- nz_cr %>%
  dplyr::filter(YearMeasured == 1) %>%
  dplyr::filter(Wave == 2018 | Wave == 2019) %>%
  dplyr::group_by(Id) %>%
  dplyr::filter(n() > 1) %>%
  dplyr::filter(n() != 0) %>%
  dplyr::ungroup(Id) %>%
  dplyr::mutate(yearS = (TSCORE - min(TSCORE)/365)) %>%
  dplyr::mutate(WSCORE = as.factor(WSCORE)) %>%
  dplyr::mutate(Covid_Timeline =
                  as.factor(ifelse(
                    TSCORE %in% 3896:3921,
                    # feb 29 - march 25th
                    "EarlyMarch",
                    ifelse(
                      TSCORE %in% 3922:3954,
                      "Lockdown",
                      #march 26- Mon 27 April 2020
                      ifelse(
                        TSCORE > 3954,
                        # after april 27th 20202
                        "PostLockdown",
                        ifelse(
                          TSCORE %in% 3842:3895,
                          # jan 6 to feb 28
                          "JanFeb",
                          ifelse(TSCORE %in% 3665:3841 &
                                   Wave == 2019,
                                 "PreCOVID",
                                 "Baseline"  # 3672 TSCORE or  20 July 2019))))))))
                          )
                        )
                      )
                    ))))%>%
  dplyr::mutate(Covid_Timeline = forcats::fct_relevel(Covid_Timeline, ord_dates_class))

Assessment

Information

Link to the NZAVS data dictionary is here

Link to questions only here

Select any one of the following three tasks

  1. Write brief report that in which you estimating demographic and ideological predictors of attitudes to government surveillance Issue.GovtSurveillance and attitudes to Issue.RegulateAI in New Zealand in 2018
Show code 
nz2018['Household.INC_s'] <- as.data.frame( scale(nz2018$Household.INC) )
nz2018['Age_s'] <- as.data.frame( scale(nz2018$Age) )
nz2018['Hours.Internet_log'] <- as.data.frame( log(nz2018$Hours.Internet + 1) )
nz2018['Issue.GovtSurveillance_s'] <- as.data.frame( log(nz2018$Issue.GovtSurveillance + 1) )
nz2018['Issue.RegulateAI_s'] <- as.data.frame( log(nz2018$Issue.RegulateAI + 1) )
nz2018['Pol.Orient_s'] <- as.data.frame( scale(nz2018$Pol.Orient) )
nz2018['Edu_s'] <- as.data.frame( scale(nz2018$Edu) )

gv <- bf(Issue.GovtSurveillance ~ Age_s + Male + Edu_s + Household.INC_s +  Pol.Orient_s)
ai <- bf(Issue.RegulateAI ~ Age_s + Male + Edu_s + Household.INC_s +  Pol.Orient_s)
fit2 <- brm(gv + ai, data = nz2018,
            file = here::here("models","workbook9_govsec_ai"))
Show code 
options(width = 120)

# Shows the prior
summary(fit2, prior = TRUE)

 Family: MV(gaussian, gaussian) 
  Links: mu = identity; sigma = identity
         mu = identity; sigma = identity 
Formula: Issue.GovtSurveillance ~ Age_s + Male + Edu_s + Household.INC_s + Pol.Orient_s 
         Issue.RegulateAI ~ Age_s + Male + Edu_s + Household.INC_s + Pol.Orient_s 
   Data: nz2018 (Number of observations: 2210) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Priors: 
Intercept_IssueGovtSurveillance ~ student_t(3, 5, 2.5)
Intercept_IssueRegulateAI ~ student_t(3, 4, 2.5)
Lrescor ~ lkj_corr_cholesky(1)
sigma_IssueGovtSurveillance ~ student_t(3, 0, 2.5)
sigma_IssueRegulateAI ~ student_t(3, 0, 2.5)

Population-Level Effects: 
                                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
IssueGovtSurveillance_Intercept           4.48      0.06     4.36     4.60 1.00     8488     2776
IssueRegulateAI_Intercept                 4.10      0.05     3.99     4.20 1.00     8363     2442
IssueGovtSurveillance_Age_s               0.08      0.04     0.01     0.16 1.00     7462     3171
IssueGovtSurveillance_MaleNot_Male        0.08      0.08    -0.07     0.23 1.00     9034     2606
IssueGovtSurveillance_Edu_s              -0.11      0.04    -0.19    -0.03 1.00     7979     3196
IssueGovtSurveillance_Household.INC_s     0.15      0.04     0.08     0.22 1.00     9364     3264
IssueGovtSurveillance_Pol.Orient_s        0.32      0.04     0.25     0.39 1.00    10387     3186
IssueRegulateAI_Age_s                     0.11      0.03     0.04     0.18 1.00     9069     3237
IssueRegulateAI_MaleNot_Male              0.53      0.07     0.40     0.67 1.00     8698     3023
IssueRegulateAI_Edu_s                    -0.08      0.03    -0.15    -0.01 1.00     8614     3102
IssueRegulateAI_Household.INC_s          -0.07      0.03    -0.13    -0.00 1.00     8357     3282
IssueRegulateAI_Pol.Orient_s              0.02      0.03    -0.05     0.09 1.00     9002     3038

Family Specific Parameters: 
                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_IssueGovtSurveillance     1.68      0.03     1.64     1.73 1.00     9831     3244
sigma_IssueRegulateAI           1.53      0.02     1.49     1.58 1.00    11920     2859

Residual Correlations: 
                                              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rescor(IssueGovtSurveillance,IssueRegulateAI)    -0.03      0.02    -0.07     0.01 1.00     9862     3196

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
brms::mcmc_plot(fit2,
               type = "areas",
               prob = .89)

Note the two outcome variables are not correlated (conditiaonl on the model):

rescor(IssueGovtSurveillance,IssueRegulateAI) -0.03

  1. Write a brief report explaining the predictors of stability and change in employment security (Emp.JobSecure) from 2018-2020 (Waves 10/11 of the the NZAVS Jitter dataset).
Show code 
ho2<- lme4::lmer(Emp.JobSecure  ~  Covid_Timeline + (1|Id), data = nzl)
plot(ggeffects::ggpredict(ho2), add.data = TRUE, jitter =.5)

$Covid_Timeline

  1. Did people gain weight during New Zealand’s Covid Lockdown?

Not a lot of evidence for it with this small sub-sample. However, this subsample might be too small.

Show code 
weight<- lme4::lmer(HLTH.BMI  ~  Covid_Timeline + (1|Id), data = nzl)
plot(ggeffects::ggpredict(weight), add.data = TRUE, dot.alpha = .1)[[1]]

Marking criteria

  1. Clarity and organisation in the descriptive component of your work.

  2. Clarity and organisation in description of your regression model.

  3. Accuracy and insight in the interpretation of your regression model.

Corrections

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Reuse

Text and figures are licensed under Creative Commons Attribution CC BY-NC-SA 4.0. Source code is available at https://go-bayes.github.io/psych-447/, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

Citation

For attribution, please cite this work as

Bulbulia (2021, May 13). Psych 447: Week 9 workbook and solutions. Retrieved from https://vuw-psych-447.netlify.app/workbooks/W_9_s/

BibTeX citation

@misc{bulbulia2021week,
  author = {Bulbulia, Joseph},
  title = {Psych 447: Week 9 workbook and solutions},
  url = {https://vuw-psych-447.netlify.app/workbooks/W_9_s/},
  year = {2021}
}