Week 10 workbook and solutions

Joseph Bulbulia https://josephbulbulia.netlify.app (Victoria University of Wellington)https://www.wgtn.ac.nz
05-23-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")
# rstan options
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores ())
theme_set(theme_classic())
Show code 
# read data

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

# to re-level 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) %>%
  dplyr::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)/365) ) %>%
  dplyr::mutate(KESSLER6sum = as.integer(KESSLER6sum))

## if you do anything with covid (warning, such a model would be tricky)

ord_dates_class <- c("Baseline",
                     "PreCOVID",
                     "JanFeb",
                     "EarlyMarch",
                     "Lockdown",
                     "PostLockdown")

nzl <- nz_cr %>%
  dplyr::filter(YearMeasured == 1) %>%
  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))%>%
  dplyr::mutate(Id = factor(Id))

Assessment

Information

Link to the NZAVS data dictionary is here

Link to questions only here

Create a distill website

  1. Create a distill website see

  2. Create a postcard for yourself [see](https://alison.rbind.io/post/2020-12-22-postcards-distill/

  3. Set up your distill webpage with your postcard on github: see

  4. Publish your postcard

Note you do not need to publish a personal website if you don’t wish to do so. We can assess you on the basis of a publishable website. Please see us for help.

Create a distill report

  1. Using the nzl dataset or another dataset of your choice, write a varying intercept/varying slope model.

Note the nzl dataset now includes 2845 individuals who have responded to at least five waves of the NZAVS. The dataset now includes relationship satisfaction Your.Personal.Relationships and employment security Emp.JobSecure

Show code 
# center religion
nzl1 <- nzl %>%
  dplyr::mutate( Relid_c = scale(Relid, scale = FALSE))%>%
  dplyr::mutate(CharityDonate = as.integer(CharityDonate))%>%
  dplyr::mutate(CharityDonate_log = log(CharityDonate + 1))

str(nzl1)

'data.frame':   22396 obs. of  83 variables:
 $ Id                         : Factor w/ 2845 levels "2","15","16",..: 1 1 1 1 1 1 1 1 2 2 ...
 $ Wave                       : Factor w/ 11 levels "2009","2010",..: 1 5 6 7 8 9 10 11 1 2 ...
 $ years                      : num  0.0411 4.1807 5.2156 6.193 7.2279 ...
 $ Age                        : num  29.7 33.7 34.8 35.7 36.8 ...
 $ Male                       : Factor w/ 2 levels "Male","Not_Male": 1 1 1 1 1 1 1 1 1 1 ...
 $ Gender                     : num  1 1 1 1 1 1 1 1 1 1 ...
 $ Edu                        : num  9 9 9 9 9 9 9 9 3 NA ...
 $ Partner                    : num  0 1 1 1 0 0 0 0 1 1 ...
 $ BornNZ                     : num  1 1 1 1 1 1 1 1 1 1 ...
 $ Employed                   : num  1 1 1 1 1 1 1 1 1 1 ...
 $ BigDoms                    : Factor w/ 5 levels "Buddhist","Christian",..: 2 4 4 4 4 4 4 4 4 4 ...
 $ TSCORE                     : num  76 1588 1966 2323 2701 ...
 $ WSCORE                     : Factor w/ 525 levels "13","14","15",..: 3 165 218 269 323 368 426 481 6 39 ...
 $ GenCohort                  : Factor w/ 5 levels "Gen Boombers: born >= 1946 & b.< 1961",..: 4 4 4 4 4 4 4 4 3 3 ...
 $ Household.INC              : num  45000 150000 110000 140000 125000 125000 130000 130000 70000 80000 ...
 $ Issue.IncomeRedistribution : num  NA NA 4 5 5 5 6 4 NA NA ...
 $ Religion.Church            : num  NA 0 0 0 0 0 0 0 NA 0 ...
 $ Religion.Believe.Cats      : num  NA 1 NA 1 1 3 3 3 NA 4 ...
 $ Your.Personal.Relationships: num  5.56 9 9 8 6 ...
 $ Census.Religion.L1         : num  2 0 0 0 0 0 0 0 0 0 ...
 $ BELONG                     : num  5.67 6.33 6.33 5.67 6 ...
 $ Religious_Group            : Factor w/ 83 levels "Anglican","Apostolic_Church_of_New_Zealand",..: 17 48 48 48 48 48 48 48 48 48 ...
 $ Relid                      : num  4 0 0 0 0 0 0 0 0 0 ...
 $ HLTH.Fatigue               : num  NA 1 1 2 1 1 1 1 NA NA ...
 $ HLTH.SleepHours            : num  NA 8 9 8 NA 8 8 9 NA NA ...
 $ HLTH.BMI                   : num  NA 23.5 24.8 24.5 23.2 ...
 $ HLTH.Weight                : num  NA 72 76 75 71 71 70 70 NA 73 ...
 $ HLTH.Height                : num  NA 1.75 1.75 1.75 1.75 1.75 1.75 1.75 NA 1.78 ...
 $ HomeOwner                  : num  NA NA NA 1 NA NA 1 NA NA NA ...
 $ Pol.Orient                 : num  NA 5 5 4 4 3 3 4 2 3 ...
 $ PATRIOT                    : num  6 7 6 5.5 6 6.5 7 6.5 5 5.5 ...
 $ Env.SatNZEnvironment       : num  8.89 8 8 8 8 ...
 $ Env.MotorwaySpend          : num  5 5 NA NA NA NA 5 5 4 NA ...
 $ Env.PubTransSubs           : num  7 6 NA NA NA NA 4 4 4 NA ...
 $ Env.ClimateChgConcern      : num  NA 1 2 2 2 5 4 3 NA NA ...
 $ LIFEMEANING                : num  NA NA NA NA NA NA 6 5.5 NA NA ...
 $ Hours.Internet             : num  NA 7 5 7 20 7 7 10 NA NA ...
 $ Issue.GovtSurveillance     : num  NA NA NA NA 5 6 7 5 NA NA ...
 $ Issue.RegulateAI           : num  NA NA NA NA NA NA 4 NA NA NA ...
 $ Hours.Exercise             : num  5 2 4 5 5 7 4 5 0 NA ...
 $ Hours.Work                 : num  40 40 40 40 40 40 45 40 40 NA ...
 $ Hours.News                 : num  NA 7 2 7 5 7 4 7 NA NA ...
 $ CONSCIENTIOUSNESS          : num  5.25 5.5 5 5.75 4.75 5.75 5.75 5.75 4 4.5 ...
 $ EXTRAVERSION               : num  5 5 4.5 5 4.5 5.5 5.25 4.5 2.75 3.5 ...
 $ AGREEABLENESS              : num  5.75 5.75 5.75 5.25 5 5.75 5.5 5.75 4.25 5 ...
 $ OPENNESS                   : num  4.5 4.25 5 5 5 5.5 5.25 5 5.5 6.25 ...
 $ Religious                  : Factor w/ 2 levels "Not_Religious",..: 2 1 1 1 1 1 1 1 1 1 ...
 $ Spiritual.Identification   : num  NA NA NA NA 5 NA 6 NA NA NA ...
 $ Believe.God                : Factor w/ 2 levels "Believe God",..: NA 1 NA 1 1 2 2 2 NA 2 ...
 $ Believe.Spirit             : Factor w/ 2 levels "Believe Spirit",..: NA 1 NA 1 1 1 1 1 NA 2 ...
 $ HoursCharity               : num  0 0 0 0 0 0 0 0 0 NA ...
 $ CharityDonate              : int  0 0 50 10 0 0 5 0 NA 20 ...
 $ Your.Future.Security       : num  7.78 8 7 8 7 ...
 $ Standard.Living            : num  8.89 9 9 9 9 ...
 $ NZ.Economic.Situation      : num  6.67 7 5 7 7 ...
 $ NZ.Social.Conditions       : num  7.78 7 7 8 8 ...
 $ NZ.Business.Conditions     : num  6.67 7 5 7 8 ...
 $ Emp.JobSecure              : num  6 7 7 6 NA 6 4 1 5 4 ...
 $ Issue.Food.GMO             : num  NA NA 4 4 5 5 5 4 NA NA ...
 $ Env.SacMade                : num  5 5 5 NA NA 6 NA NA 5 6 ...
 $ KESSLER6sum                : int  NA 3 3 5 2 5 2 3 NA 4 ...
 $ YearMeasured               : num  1 1 1 1 1 1 1 1 1 1 ...
 $ Warm.Muslims               : num  NA 2 1 3 3 4 2 4 NA NA ...
 $ Warm.Immigrants            : num  4 3 4 4 4 4 4 4 4 4 ...
 $ FeelHopeless               : Factor w/ 5 levels "None Of The Time",..: NA 2 2 2 1 1 1 2 NA 1 ...
 $ FeelDepressed              : Factor w/ 5 levels "None Of The Time",..: NA 1 1 1 1 1 1 1 NA 1 ...
 $ FeelRestless               : Factor w/ 5 levels "None Of The Time",..: NA 2 2 2 2 3 2 2 NA 3 ...
 $ EverythingIsEffort         : Factor w/ 5 levels "None Of The Time",..: NA 1 1 3 1 2 1 1 NA 1 ...
 $ FeelWorthless              : Factor w/ 5 levels "None Of The Time",..: NA 1 1 1 1 1 1 1 NA 1 ...
 $ FeelNervous                : Factor w/ 5 levels "None Of The Time",..: NA 2 2 2 2 3 2 2 NA 3 ...
 $ male_id                    : Factor w/ 2 levels "Male","Not_Male": 1 1 1 1 1 1 1 1 1 1 ...
 $ date                       : Date, format: "2009-09-14" ...
 $ FeelWorthless_int          : int  NA 1 1 1 1 1 1 1 NA 1 ...
 $ FeelNervous_int            : int  NA 2 2 2 2 3 2 2 NA 3 ...
 $ FeelHopeless_int           : int  NA 2 2 2 1 1 1 2 NA 1 ...
 $ EverythingIsEffort_int     : int  NA 1 1 3 1 2 1 1 NA 1 ...
 $ FeelRestless_int           : int  NA 2 2 2 2 3 2 2 NA 3 ...
 $ FeelDepressed_int          : int  NA 1 1 1 1 1 1 1 NA 1 ...
 $ HLTH.Fatigue_int           : int  NA 2 2 3 2 2 2 2 NA NA ...
 $ yearS                      : num  75.8 1587.8 1965.8 2322.8 2700.8 ...
 $ Covid_Timeline             : Factor w/ 6 levels "Baseline","PreCOVID",..: 1 1 1 1 1 1 1 2 1 1 ...
 $ Relid_c                    : num [1:22396, 1] 2.13 -1.87 -1.87 -1.87 -1.87 ...
  ..- attr(*, "scaled:center")= num 1.87
 $ CharityDonate_log          : num  0 0 3.93 2.4 0 ...
Show code 

# demean
nzl1 <- cbind(nzl1, 
           parameters::demean(nzl1, 
                  c("Relid"),
                  group = "Id")
           )

Model:

Show code 
nzl1 <- nzl1 %>%
  dplyr::mutate(Relid_between_c = scale(Relid_between, center = TRUE, scale = FALSE))  
# test dataset
m1 <- brms::brm(CharityDonate ~ yearS + Relid_c  + (1 + Relid_c|Id),  
                family = "negbinomial",
                file = here::here("models", "workbook_10_ex1"),
                data = nzl1)

summary(m1)

Note poor mixing. This is a sign of a problem:

Show code 
# results
parameters::model_parameters(
  m1,
  dispersion = FALSE,
  test = "pd",
  priors = TRUE,
  effects = "all"
)%>% print_md()
Show code 
# graph
p_m1 <- plot(brms::conditional_effects(m1,
                            "Relid_c",
                            spaghetti = TRUE,
                            nsamples = 100)
             )[[1]]
p_m1

Attempt model with demeaned data

Show code 
m2 <-
  brms::brm(
    CharityDonate ~ yearS + Relid_within + Relid_between_c  + 
      (1 + Relid_within | Id),
    family = "negbinomial",
    file = here::here("models", "workbook_10_ex2"),
    data = nzl1
  )

We can see that people are rounding in their estimates:

Show code 
summary(m2, prior = TRUE)
brms::pp_check(m2) + scale_x_continuous(limits = c(0,10000))
Show code 
# graph
p_m2a <- plot(brms::conditional_effects(m2,
                            c("Relid_within"),
                            spaghetti = TRUE,
                            re_formula = NULL,
                            nsamples = 100)
             )
p_m2a
Show code 
p_m2b <- plot(brms::conditional_effects(m2,
                            c("Relid_between_c"),
                            spaghetti = TRUE,
                            re_formula = NULL, # all group-level effects
                            nsamples = 100)
             )
p_m2b

Log transform outcomes

Show code 
# Response distributions in brms
#https://cran.r-project.org/web/packages/brms/vignettes/brms_families.html

set.seed(123)
nm <- sample(nzl1$Id, size = 1000)
length(unique(nzl1$Id))

[1] 2845
Show code 

nzl2 <- nzl1 %>%
 filter(Id %in% nm) %>%
  droplevels()

# 825
length(unique(nzl2$Id))

[1] 825
Show code 
  
# What is this? 
hist(nzl1$CharityDonate_log)
Show code 
# sum(nzl1$CharityDonate_log == 0, na.rm=TRUE)/nrow(nzl1) # 14% zero
m3 <-
  brms::brm(
    CharityDonate_log ~  yearS + Relid_within + Relid_between_c  +
      (1 + Relid_within | Id),
    family = "hurdle_lognormal",
    file = here::here("models", "workbook_10_ex3"),
    data = nzl2
  )

Mixing is not great, model should be run for longer, however

Show code 
summary(m3, prior = TRUE)
brms::pp_check(m3)# + scale_x_continuous(limits = c(0,10000))

Within person effects are not, on this model, sample, clear.

Show code 
# graph
p_m3a <- plot(brms::conditional_effects(m3,
                            c("Relid_within"),
                            spaghetti = TRUE,
                            re_formula = NULL,
                            nsamples = 100)
             )
Show code 
p_m3a

$Relid_within

The religion charity effect is evident

Show code 
p_m3b <- plot(brms::conditional_effects(m3,
                            c("Relid_between_c"),
                            spaghetti = TRUE,
                            re_formula = NULL, # all group-level effects
                            nsamples = 100)
             )
Show code 
p_m3b

$Relid_between_c

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.

Reports I’ve been doing FYI

https://go-bayes.github.io/reports/

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 23). Psych 447: Week 10 workbook and solutions. Retrieved from https://vuw-psych-447.netlify.app/workbooks/W_10_s/

BibTeX citation

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