Import the jittered NZAVS dataset, selecting the 2018 Wave
### Libraries
library("tidyverse")
library("patchwork")
library("lubridate")
library("kableExtra")
library("gtsummary")
library("lubridate")
#devtools::install_github("data-edu/tidyLPA")
if (!require(tidyLPA)) {
install.packages("tidyLPA")
}
# 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 <- 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)) %>%
mutate(FeelHopeless = forcats::fct_relevel(FeelHopeless, f)) %>%
mutate(FeelDepressed = forcats::fct_relevel(FeelDepressed, f)) %>%
mutate(FeelRestless = forcats::fct_relevel(FeelRestless, f)) %>%
mutate(EverythingIsEffort = forcats::fct_relevel(EverythingIsEffort, f)) %>%
mutate(FeelWorthless = forcats::fct_relevel(FeelWorthless, f)) %>%
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::filter(Wave == 2018)
dplyr::glimpse(nz)
Rows: 2,918
Columns: 82
$ Id <dbl> 1, 4, 5, 6, 7, 9, 12, …
$ Wave <fct> 2018, 2018, 2018, 2018…
$ years <dbl> 9.921971, 9.900068, 9.…
$ Age <dbl> 43.85489, 47.49350, 33…
$ Male <fct> Not_Male, Male, Not_Ma…
$ Gender <dbl> 0, 1, 0, 1, 0, 0, 0, 1…
$ w.GendAgeEthnic <dbl> 0.4022481, 4.1207264, …
$ Edu <dbl> 7, 7, 7, 7, 7, 1, 4, 9…
$ Partner <dbl> 1, 1, 1, 1, 0, 0, 1, 1…
$ BornNZ <dbl> 1, 1, 1, 1, 1, 1, 1, 1…
$ Employed <dbl> 1, 1, 1, 1, 1, 0, 0, 1…
$ BigDoms <fct> Not_Rel, Not_Rel, Not_…
$ TSCORE <dbl> 3685, 3677, 3687, 3399…
$ GenCohort <fct> GenX: born >=1961 & b.…
$ Hours.Exercise <dbl> 2, 10, 1, 3, 10, 4, 0,…
$ Hours.Work <dbl> 16.0, 65.0, 0.0, 35.0,…
$ Hours.News <dbl> 2.0, 4.0, 0.5, 5.0, 1.…
$ Hours.Internet <dbl> 7, 15, 10, 20, 7, 0, 5…
$ Hours.SocialMedia <dbl> 4, 5, 10, 2, 2, 0, 2, …
$ HoursCharity <dbl> 6, 2, 3, 0, 0, 0, 0, 0…
$ Hours.CompGames <dbl> 0, 0, 5, 2, 0, 0, 2, 0…
$ Hours.Family <dbl> 7.00, 1.00, 10.00, 0.0…
$ Hours.Friends <dbl> 4, 3, 5, 2, 10, 0, 8, …
$ Hours.Community <dbl> 0, 0, 2, 0, 0, 0, 4, 0…
$ Hours.Religious <dbl> 0, 5, 0, 0, 0, 0, 0, 0…
$ CharityDonate <dbl> 200, 100, 200, 0, 50, …
$ Family.Money <dbl> 0, 0, 0, 0, 200, 0, 0,…
$ Friends.Money <dbl> 0, 0, 0, 0, 0, 0, 0, 0…
$ Community.Money <dbl> 0, 0, 0, 0, 0, 0, 0, 0…
$ Family.Time <dbl> 0, 0, 12, 0, 0, 0, 0, …
$ Friends.Time <dbl> 0, 0, 1, 0, 0, 0, 0, 0…
$ Community.Time <dbl> 0, 0, 0, 0, 0, 0, 0, 0…
$ Household.INC <dbl> 200000, 40000, 140000,…
$ Issue.IncomeRedistribution <dbl> 6, 4, 4, 7, 7, 4, 3, 5…
$ Religion.Church <dbl> 0, 0, 0, 0, 0, 0, 0, 0…
$ Religion.Believe.Cats <dbl> 4, 1, 4, 1, 1, 1, 3, 3…
$ Relid <dbl> 0, 0, 0, 5, 0, 4, 0, 0…
$ HLTH.Fatigue <dbl> 1, 2, 2, 1, 1, 1, 3, 1…
$ HLTH.SleepHours <dbl> 7.0, 4.0, 6.0, 7.0, 7.…
$ HLTH.BMI <dbl> 19.83471, 13.14828, 26…
$ HLTH.Weight <dbl> 54.0, 45.0, 70.0, 70.0…
$ HLTH.Height <dbl> 1.65, 1.85, 1.63, 1.80…
$ HomeOwner <dbl> 1, 0, 1, 0, 1, 0, 1, 1…
$ Pol.Orient <dbl> 5, 3, 4, 1, 2, 4, 7, 2…
$ PATRIOT <dbl> 6.0, 7.0, 7.0, 7.0, 6.…
$ Env.SatNZEnvironment <dbl> 7, 7, 9, 3, 3, 3, 6, 6…
$ Env.MotorwaySpend <dbl> 3, 5, 4, 3, 5, 7, 5, 6…
$ Env.PubTransSubs <dbl> 7, 5, 6, 7, 6, 4, 7, 5…
$ Env.ClimateChgConcern <dbl> 6, 7, 5, 7, 7, NA, 7, …
$ LIFEMEANING <dbl> 6.0, 4.5, 5.0, 7.0, 7.…
$ Issue.GovtSurveillance <dbl> 5, 3, 4, 7, 6, 7, 6, 4…
$ Issue.RegulateAI <dbl> 6, 4, 5, 4, 2, 5, 4, 6…
$ CONSCIENTIOUSNESS <dbl> 5.250000, 5.500000, 2.…
$ EXTRAVERSION <dbl> 4.00, 4.00, 4.75, 2.25…
$ AGREEABLENESS <dbl> 6.75, 6.00, 6.00, 4.25…
$ OPENNESS <dbl> 2.00, 4.25, 4.25, 6.25…
$ Religious <fct> Not_Religious, Not_Rel…
$ Spiritual.Identification <dbl> 2, 5, 2, 7, 7, NA, 4, …
$ Believe.God <fct> Not Believe God, Belie…
$ Believe.Spirit <fct> Not Believe Spirit, Be…
$ Your.Personal.Relationships <dbl> 10, 2, 7, 10, 8, 10, 2…
$ Your.Future.Security <dbl> 10, 6, 8, 7, 8, 8, 2, …
$ Standard.Living <dbl> 10, 6, 9, 10, 8, 7, 10…
$ NZ.Economic.Situation <dbl> 8, 6, 4, 10, 3, 2, 3, …
$ NZ.Social.Conditions <dbl> 7, 6, 3, 4, 3, 3, 0, 7…
$ NZ.Business.Conditions <dbl> 8, 6, 5, 10, 3, 6, 2, …
$ Emp.JobSecure <dbl> 7, 6, 5, 7, 6, NA, NA,…
$ Emp.JobSat <dbl> 6, 6, 5, 7, 6, NA, NA,…
$ Emp.JobValued <dbl> 7, 5, 5, 7, 6, NA, NA,…
$ Issue.Food.GMO <dbl> 6, 5, 6, 7, 7, 7, 4, 5…
$ Env.SacMade <dbl> NA, NA, NA, NA, NA, NA…
$ KESSLER6sum <dbl> 4, 7, 10, 2, 3, 2, 8, …
$ POLICE.ENGAGE <dbl> 7.0, 4.5, 7.0, 6.5, 6.…
$ POLICE.TRUST <dbl> 5.666667, 4.333333, 5.…
$ FeelHopeless <fct> A Little Of The Time, …
$ FeelDepressed <fct> None Of The Time, None…
$ FeelRestless <fct> A Little Of The Time, …
$ EverythingIsEffort <fct> A Little Of The Time, …
$ FeelWorthless <fct> None Of The Time, None…
$ FeelNervous <fct> A Little Of The Time, …
$ male_id <fct> Not_Male, Male, Not_Ma…
$ date <date> 2019-08-02, 2019-07-2…Assessment
Information
Link to the NZAVS data dictionary is here
Link to questions only here
Task
- Write brief report that predicts belief in spirit or a life-force (
Believe.Spirit) from no more than five covariates. Explain your model and results.
Approach to the solutions
Note that you’ll need to relevel your key factor.
Also, for a Poisson outcome, the indicator needs to be an integer.
# Relevel factor
library(ggplot2)
nz1 <- nz %>%
filter(Wave == 2018) %>%
mutate(Believe.Spirit = forcats::fct_relevel(Believe.Spirit, c("Not Believe Spirit",'Believe Spirit'))) %>%
mutate(Believe.God = forcats::fct_relevel(Believe.God, c("Not Believe God",'Believe God'))) %>%
mutate(CharityDonate = as.integer(CharityDonate))
Because religion is often a non-linear predictor, I’d attempt a spline model:
library("splines")
m1<-glm(Believe.Spirit ~ bs(Relid), data = nz1, family = binomial())
p <- plot(
ggeffects::ggpredict(m1)[[1]],
)
p + scale_y_continuous(limits = c(0,1)) +
xlab("Religious Identification 1-7") +
theme_classic()
Tapering of credulity among the highly religious. Why?
- Write report that predicts charitable donations (
CharityDonate) using no more than five co-variates. Explain your model and results.
Inspect data:
hist(nz1$CharityDonate,breaks=1000)
For this model, I would again use a spline:
library(splines)
m2<-glm(CharityDonate ~ bs(Relid), data = nz1, family = gaussian())
p2 <- plot(
ggeffects::ggpredict(m2, effects = "Relid[all]")[[1]],
add.data = TRUE,
dot.alpha = .2
)
p2 +
xlab("Religious Identification 1-7") +
#scale_y_log10() +
theme_classic()
Consider following graph tells a better story (why, why not?)
p3 <- plot(
ggeffects::ggpredict(m2)[[1]],
)
p3 +
xlab("Religious Identification 1-7") +
theme_classic()
Compare different families. Compare a model without splines. Polish the graphs… tell the story.
Marking criteria
Clarity and organisation in the descriptive component of your work.
Clarity and organisation in description of your regression model.
Accuracy and insight in the interpretation of your regression model.