Financial satisfaction and Religion

Introduction

Financial satisfaction reflects how individuals perceive their financial situation and their capacity to meet personal goals. While objective income is a clear determinant, subjective and psychosocial factors are equally consequential for well-being. Among these, religiosity has emerged as a significant yet debated predictor.

The Complex Link Between Faith and Financial Well-being

The relationship between religiosity and financial satisfaction is complex, context-dependent, and often contradictory (Cwynar et al., 2024; Kose & Cinar, 2024). Higher religiosity frequently correlates with greater financial well-being, yet no consensus exists on a single directional effect. While Sarofim et al. (2020) documented a positive influence of religion on financial well-being, Migheli (2018) found that specific religious groups report lower satisfaction than secular peers. Cross-national surveys further complicate the picture, revealing a nonlinear, U- or J-shaped pattern in which modest religiousness may slightly depress financial satisfaction, whereas high religiosity yields a strongly positive effect (Kose & Cinar, 2024).

The Dual Mechanisms of Religiosity: Belief and Attendance

Unpacking this complexity requires operationalizing religiosity through two core dimensions: intrinsic belief in God and organizational attendance (Kose & Cinar, 2024). Each offers a distinct pathway to financial satisfaction. Intrinsic faith provides meaning and hope, functioning as a psychological coping mechanism that buffers financial anxiety (Berkessel et al., 2021). Many faiths also promote modesty and discourage materialism, reducing pressure for conspicuous spending and thereby increasing reported satisfaction at a given income level (Friedman & Lynch, 2023; Kalbasi & Amani, 2022). Religious attendance, in turn, fosters robust social networks that generate social capital, interpersonal trust, and mutual aid — resources that can directly improve financial behaviors such as saving (Lim et al., 2010; Kortt et al., 2015; Hong et al., 2024).

The Present Study

We argue that inconsistencies in the literature stem from two sources: the conflation of distinct religiosity dimensions, and the neglect of country-level contextual variation. This study addresses both. Regarding dimensionality, belief in God and religious attendance likely operate through distinct psychological and social mechanisms whose effects on financial satisfaction need not be equivalent or even directionally consistent. Belief primarily functions through meaning-making and cognitive reappraisal of material circumstances, whereas attendance provides access to tangible social resources — community networks, reciprocal support, and collective identity.

Accordingly:

H1: Belief in God and religious attendance will show differential effects on financial satisfaction, with religious attendance exhibiting a stronger positive association due to its social resource dimension.

Regarding contextual variation, we draw on two complementary theoretical frameworks. Existential Security Theory (Norris & Inglehart, 2004) holds that religion’s psychological and material functions are contingent on societal vulnerability. In advanced welfare states, secular institutions provide existential security that substitutes for the coping resources religion traditionally offers, attenuating its buffering effect. In developing nations with weaker safety nets, religion remains a crucial resource for managing risk and offsetting material precarity.

H2: The positive association between religiosity and financial satisfaction will be stronger in countries with lower economic security, and weaker or absent in advanced welfare states, reflecting the substitution of secular institutions for religious coping resources.

The religiosity-as-social-value hypothesis (Gebauer et al., 2017) offers a complementary account. Where religion is culturally normative, religious participation confers social approval, status, and inclusion — amplifying the individual benefits of religiosity not through any intrinsic property of faith, but through conformity to locally valued identities.

H3: Countries with higher average belief in God will exhibit a stronger positive association between religion and financial satisfaction.

H4: This contextual amplification by national religiosity will be weaker or absent for religious attendance, given its more behaviorally constrained nature compared to belief.

Exploratory Analysis

In this section we will mainly select and describe the data needed to conduct this current study. We will work with the world values survey (wvs) which will allow us to grasp the difference between countries which is the main goal for this study.

Variables selection

wvs <- read.csv( '/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Multilevel Analysis/Final Project/WVS.csv', header = T) 

wvs <- wvs %>% select(country = B_COUNTRY_ALPHA, finsat = Q50, Incme_lvl = Q288, age = Q262, sex = Q260, edu = Q275, belief_god = Q164, religious_at = Q171) %>% filter(!country == "IRQ")


hdi <- readxl::read_excel("~/Documents/Old Pc/KUL QASS/Multilevel Analysis/Final Project/HDR25_Statistical_Annex_HDI_Table.xlsx", skip = 5) %>% filter(!is.na(`HDI rank`))


wvs$country_name <- countrycode(sourcevar = wvs$country, 
                                origin = "iso3c", 
                                destination = "country.name",
                                custom_match = c("NIR" = "Northern Ireland"))


hdi <- hdi %>% select(country = Country, Value)


hdi$iso3c <- countrycode(sourcevar = hdi$country, 
                         origin = "country.name", 
                         destination = "iso3c")

wvs$iso3c <- countrycode(sourcevar = wvs$country_name, # or wvs$country
                         origin = "country.name", 
                         destination = "iso3c",
                         custom_match = c("Northern Ireland" = "NIR"))


wvs <- wvs %>%
  mutate(iso3c = ifelse(country_name == "Northern Ireland", "GBR", iso3c))

df <- left_join(wvs, hdi, by = "iso3c")

df <- df %>% mutate(hdi = Value) %>% select(-Value)

df <- df %>% group_by(country_name) %>% mutate(mean_belief_god = mean(belief_god, na.rm = TRUE)
                                                   ) %>%  ungroup()

We combine data from the World Values Survey with data from the Human Development Index to produce the complete DataFrame for this study.

Variables Transformation

#table(f_df$finsat)
df$religious_at <-  8 - ifelse(df$religious_at < 0, NA, df$religious_at) 

df$finsat <- ifelse(df$finsat < 0, NA, df$finsat)

df$Incme_lvl <- ifelse(df$Incme_lvl < 0, NA, df$Incme_lvl)


df$age <- ifelse(df$age < 0, NA, df$age)


df$sex <- ifelse(df$sex < 0, NA, df$sex)
df$sex <- as.factor(df$sex)
df$sex <- ifelse(df$sex  == "1", "Male", "Female")
df$sex <-  as.factor(df$sex)


df$edu <- ifelse(df$edu < 0, NA, df$edu)


df$belief_god <- ifelse(df$belief_god < 0, NA, df$belief_god)



df$hdi <- round(as.numeric(df$hdi), 2)
df <- df %>%
  mutate(hdi = case_when(
    country_name == "Taiwan" ~ 0.93,
    country_name == "Puerto Rico" ~ 0.87,
    country_name == "Macao SAR China" ~ 0.93,
    TRUE ~ hdi # If it's not one of those three, keep the existing HDI value
  ))

Below are our variables after transformation

Label <- c("Fin_sat","country", "religious_at", "belief_god","Incme_lvl", "age", "edu", "sex",  "hdi")

Meaning <- c("Financial satisfaction","Country Name", "Religious Services Attendance", "Belief in the Importance of God", "Income Level", "Age", "Education Level","Sex", "Human Development Index")

Level_Of_Measurement <- c("(Quasi-)Interval", "Ordinal", "(Quasi-)Interval", "(Quasi-)Interval", "(Quasi-)Interval", "Ratio", "Categorical", "Categorical", "Ratio")

R_Data_Type <- c(class(df$finsat), class(df$country_name), class(df$religious_at), class(df$belief_god), class(df$Incme_lvl), class(df$age), class(df$edu), class(df$sex), class(df$hdi))

df_val <- data.frame(Label, Meaning, Level_Of_Measurement, R_Data_Type, stringsAsFactors = FALSE)  

kable(df_val) %>% 
  kable_styling(bootstrap_options=c("bordered", "responsive","striped"), full_width = FALSE)

Label	Meaning	Level_Of_Measurement	R_Data_Type
Fin_sat	Financial satisfaction	(Quasi-)Interval	integer
country	Country Name	Ordinal	character
religious_at	Religious Services Attendance	(Quasi-)Interval	numeric
belief_god	Belief in the Importance of God	(Quasi-)Interval	integer
Incme_lvl	Income Level	(Quasi-)Interval	integer
age	Age	Ratio	integer
edu	Education Level	Categorical	integer
sex	Sex	Categorical	factor
hdi	Human Development Index	Ratio	numeric

Below are the tables summarising our variables.

df <- df %>% select(-country.x, -iso3c, -country.y)
skim(df)

Data summary
Name	df
Number of rows	93078
Number of columns	10
_______________________
Column type frequency:
character	1
factor	1
numeric	8
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
country_name	0	1	4	19	0	63	0

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
sex	95	1	FALSE	2	Fem: 49188, Mal: 43795

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
finsat	578	0.99	6.22	2.42	1.00	5.00	6.00	8.00	10.00	▂▃▇▇▅
Incme_lvl	2928	0.97	4.90	2.08	1.00	4.00	5.00	6.00	10.00	▃▅▇▃▁
age	510	0.99	43.50	16.60	16.00	30.00	42.00	56.00	103.00	▇▇▆▂▁
edu	1009	0.99	3.56	2.02	0.00	2.00	3.00	6.00	8.00	▃▇▂▅▂
belief_god	1004	0.99	7.33	3.25	1.00	5.00	9.00	10.00	10.00	▂▁▂▂▇
religious_at	1062	0.99	3.77	2.20	1.00	1.00	4.00	6.00	7.00	▇▁▃▂▆
hdi	0	1.00	0.81	0.11	0.50	0.73	0.80	0.93	0.96	▁▁▆▆▇
mean_belief_god	0	1.00	7.23	2.21	2.79	4.93	7.85	9.23	9.92	▂▃▃▂▇

We have in total 93078 observations across 63 countries.

In the next section we will explore the distribution of each variables, to have a overall understanding of our sample.

Evaluation and Description of Distribution

In this section we will closely examine the distribution of each variable to get a first insight into our data, and perceive the flaws which might, in further sections, impede the quality of our model.

Financial satisfaction Distribution

densityplot(~df$finsat | df$country_name) # distribution of age conditional on country = by country

Results: In most countries, people are more or less satisfied with their financial situation, which can be seen from the high prevalence of the left-tail distribution. However, some countries, such as Brazil and Colombia, have a right-tail distribution, suggesting that most people are generally less satisfied with their financial situation in these countries.

Religious Attendance

densityplot(~df$religious_at | df$country_name)

Results: Religious services attendance highly varies between countries. For instance, in Canada most of the respondents less often or do not attend at all religious services while in Indonesia, the majority of respondent regularly attend religious services.

Importance of God

densityplot(~df$belief_god| df$country_name)

Results: From one country to another, belief in the importance of God is evenly distributed across the scores. But for some countries, such as Egypt and Bolivia, there is a high density of respondents who believe extremely strongly in the importance of God in their lives.

Missing Values and Imputation

Running a model with missing values can prevent us from obtaining a complete picture of our research. Indeed, many predictive models simply remove observations containing missing values before running them. This considerably reduces the number of observations and leads to potential flaws that may influence the significance of our study. It is therefore important to deal with missing values before proceeding with our study.

First we look at variables containing missing values.

unique_countries <- unique(df$country_name)

df$country_name <- as.factor(df$country_name)
country_dfs <- list()

for (country in unique_countries) {
  country_df <- df[df$country_name == country, ]
  country_dfs[[country]] <- as.data.frame(country_df)
}

sapply(df, function(x) sum(is.na(x)))

##          finsat       Incme_lvl             age             sex             edu 
##             578            2928             510              95            1009 
##      belief_god    religious_at    country_name             hdi mean_belief_god 
##            1004            1062               0               0               0

As we can see, all our first-level variables contain missing values, with income level having the highest prevalence (2928 missing values), followed by religious attendance (1062 missing values). Income related questions are knwon to be the most sensitive questions in a survey. sensitive questions.

Next we look at the pattern of missingness across countries to determine whether those missing are MAR, MCAR, or MNAR.

Pattern of missingness

mdf <- list()

for (country in  unique_countries) {
  
  mdf[[country]] <- missing_data.frame (as.data.frame(country_dfs[[country]]) )
  
  mdf[[country]] <- change (mdf[[country]] , y = c("age", "finsat", "Incme_lvl", "belief_god", "life_sat", "religious_at", "hdi"), what = "type", to = "positive-continuous")
}

(We had a look at all the missing images per country). For most of the the countries, the NAs have no particular pattern that may suggest that they are MNAR. But for some countries, such as New Zealand there seems to be a certain pattern of missing data. We should probably take a closer look at them.

New Zealand missing patterns

unique_countries[39]

## [1] "New Zealand"

image(mdf[[39]])

We can see here that there is an important pattern of missingness around the level of income. we should analyse that particular variables across other variables explaining variables.

 country_dfs[[39]] %>% mosaicMiss(highlight = "Incme_lvl", plotvars = c("edu"))

Analysis of the data distribution reveals a systematic pattern where the proportion of missing income records decreases sequentially as participants’ higher education levels increase. This suggests that the income data are Missing at Random (MAR) conditional on education, necessitating the inclusion of the education variable in subsequent imputation procedures to ensure unbiased estimates.

Imputation

For this project, we decided to impute missing data using the MissForests method. As we saw earlier, patterns of missing data vary greatly from country to country. It would therefore be more appropriate to impute missing values separately for each country. Below is the commented code for our imputation. As the country imputation process takes too long, we decided to run it once and save the imputed database in a csv file which we reloaded for further processing.

registerDoParallel(cores=4)

imputed_missForests <- list()

for  (country in unique_countries){
    imputed_missForests[[country]] <- missForest(country_dfs[[country]], verbose = TRUE, parallelize = "forests")
  }
 imputed_missForests[[1]]$OOBerror[1]

missForest_imputed <- data.frame(
   original = country_dfs[[1]]$finsat,
  imputed_missForest = imputed_missForests[[1]]$ximp$finsat
  )

plot_data <- missForest_imputed %>%
pivot_longer(cols = c(original, imputed_missForest),
            names_to = "Data_Type",
             values_to = "Financial_Satisfaction")

# 2. Create an overlaid density plot
ggplot(plot_data, aes(x = Financial_Satisfaction, fill = Data_Type)) +
 geom_density(alpha = 0.5) + # alpha controls the transparency
theme_minimal() +
labs(title = "Distribution of Financial Satisfaction",
subtitle = "Original Data vs. missForest Imputation",
x = "Financial Satisfaction Score",
 y = "Density",
fill = "Legend") +
scale_fill_manual(values = c("original" = "#1f77b4", "imputed_missForest" = "#ff7f0e"),
              labels = c("Imputed (missForest)", "Original (with NAs)"))

clean_df <- data.frame()

for  (country in unique_countries){
 clean_df <- rbind(clean_df, as.data.frame(imputed_missForests[[country]]$ximp))
 }
sapply(imputed_missForests, function(x) ncol(x$ximp))

data <- write.csv(clean_df, '/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Multilevel Analysis/Final Project/data5.csv')

As we can see the miss forest imputated financial satisfaction data follows the same distribution than the original one.

Imputed Dataset

Next, we imported our new imputed dataset.

data <- read.csv('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Multilevel Analysis/Final Project/data5.csv', header = T) %>% dplyr::select(-X)


data <- data %>% group_by(country_name) %>% mutate(mean_belief_god = mean(belief_god, na.rm = TRUE)
                                                   ) %>%  ungroup()


data$belief_god_centred <- center(data$belief_god, type ="CWC", cluster = data$country_name)


#head(data)



data <-  data %>% dplyr::select(country = country_name, finsat, belief_god, religious_attendance = religious_at, income_level=Incme_lvl, age, sex, education_level = edu,hdi,mean_belief_god , belief_god_centred)

sapply(data, function (x) sum(is.na(data)))

##              country               finsat           belief_god 
##                    0                    0                    0 
## religious_attendance         income_level                  age 
##                    0                    0                    0 
##                  sex      education_level                  hdi 
##                    0                    0                    0 
##      mean_belief_god   belief_god_centred 
##                    0                    0

#isf_df$avg_life_sat


# #head(f_df)
# 
# f_df <- f_df %>% select(country = country_names, finsat, Incme_lvl, householdsize , age, sex,edu,avg_life_sat, inflation.rate, Local.Purchasing.Power.Index, ginindex)

As can be seen, the newly imputed dataset has no missing values. It now contains the first and second level variables required for the current study. In the next section, we performed a bivariate analysis between our IVs and our DV to test whether the linearity assumption is fulfilled to run our model.

skim(data) %>% select(-n_missing, -complete_rate) %>% 
  mutate(across(where(is.numeric), ~ round(.x, 2)))

Data summary
Name	data
Number of rows	93078
Number of columns	11
_______________________
Column type frequency:
character	2
numeric	9
________________________
Group variables	None

Variable type: character

skim_variable	min	max	empty	n_unique	whitespace
country	4	19	0	63	0
sex	4	6	0	2	0

Variable type: numeric

skim_variable	mean	sd	p0	p25	p50	p75	p100	hist
finsat	6.22	2.41	1.00	5.00	6.00	8.00	10.00	▂▃▇▇▅
belief_god	7.31	3.25	1.00	5.00	9.00	10.00	10.00	▂▁▂▂▇
religious_attendance	3.76	2.19	1.00	1.00	4.00	6.00	7.00	▇▁▃▂▆
income_level	4.90	2.06	1.00	4.00	5.00	6.00	10.00	▃▅▇▃▁
age	43.51	16.57	16.00	30.00	42.00	56.00	103.00	▇▇▆▂▁
education_level	3.57	2.01	0.00	2.00	3.00	5.16	8.00	▃▇▂▅▂
hdi	0.81	0.11	0.50	0.73	0.80	0.93	0.96	▁▁▆▆▇
mean_belief_god	7.31	2.16	2.80	5.04	7.92	9.29	9.96	▂▃▃▂▇
belief_god_centred	0.00	2.43	-8.88	-1.28	0.37	1.02	7.20	▁▂▇▅▁

table(data$sex)

## 
## Female   Male 
##  49243  43835

Bivariates Tests

Financial satisfaction vs Numeric IVs

First we look at the relationship between our DV and the selected continuous IVs.

data <-  data %>% 
  mutate(across(where(is.numeric), ~scale(.)))

sjPlot::tab_corr(data %>% dplyr::select(-country, -sex), corr.method = "pearson")

	finsat	belief_god	religious_attendance	income_level	age	education_level	hdi	mean_belief_god	belief_god_centred
finsat		-0.028***	-0.029***	0.335***	0.018***	0.106***	0.079***	-0.081***	0.035***
belief_god	-0.028***		0.542***	-0.067***	-0.088***	-0.185***	-0.474***	0.665***	0.747***
religious_attendance	-0.029***	0.542***		-0.036***	-0.053***	-0.155***	-0.387***	0.478***	0.300***
income_level	0.335***	-0.067***	-0.036***		-0.098***	0.286***	0.096***	-0.072***	-0.026***
age	0.018***	-0.088***	-0.053***	-0.098***		-0.123***	0.279***	-0.234***	0.091***
education_level	0.106***	-0.185***	-0.155***	0.286***	-0.123***		0.283***	-0.219***	-0.052***
hdi	0.079***	-0.474***	-0.387***	0.096***	0.279***	0.283***		-0.713***	0.000
mean_belief_god	-0.081***	0.665***	0.478***	-0.072***	-0.234***	-0.219***	-0.713***		0.000
belief_god_centred	0.035***	0.747***	0.300***	-0.026***	0.091***	-0.052***	0.000	0.000
Computed correlation used pearson-method with listwise-deletion.

Financial satisfaction vs Categorical IVs

Financial Satisfaction vs Sex

t.test(data$finsat ~ data$sex)

## 
##  Welch Two Sample t-test
## 
## data:  data$finsat by data$sex
## t = -2.7963, df = 92086, p-value = 0.005171
## alternative hypothesis: true difference in means between group Female and group Male is not equal to 0
## 95 percent confidence interval:
##  -0.03120827 -0.00548731
## sample estimates:
## mean in group Female   mean in group Male 
##         -0.008640875          0.009706915

plotmeans(finsat ~ sex, data = data, 
          xlab = "Sex", ylab = "Household financial satisfaction",
          main="Mean Plot with 95% CI")

Multilevel regression

fit <- lm(finsat ~ 1, data = data)
nullmodel <- lmer(finsat ~ 1 +  (1 | country), data = data, REML = FALSE) 
summary(nullmodel)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: finsat ~ 1 + (1 | country)
##    Data: data
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##  255430.4  255458.7 -127712.2  255424.4     93075 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.81243 -0.63822  0.04849  0.71806  2.75573 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  country  (Intercept) 0.1049   0.3240  
##  Residual             0.9074   0.9526  
## Number of obs: 93078, groups:  country, 63
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -0.01921    0.04095  -0.469

2*(logLik(nullmodel) - logLik(fit))

## 'log Lik.' 8718.533 (df=3)

icc(nullmodel)

## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.104
##   Unadjusted ICC: 0.104

Based on the ICC of the null model, ten percent of the overall variance in our dataset can be explained by the second-level variable, suggesting that it makes sense for this study to apply hierarchical regression with country as a random intercept. Furthermore, the random effect suggests that there is variance in financial position within and between countries.

Forward Selection Models

m1x <- lmer(finsat ~  religious_attendance + (1 | country), data = data, REML = FALSE)
m2x <- lmer(finsat ~  religious_attendance + belief_god_centred +  (1 | country), data = data, REML = FALSE)

anova(m1x, m2x, test = "Chisq")

## Data: data
## Models:
## m1x: finsat ~ religious_attendance + (1 | country)
## m2x: finsat ~ religious_attendance + belief_god_centred + (1 | country)
##     npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
## m1x    4 255308 255346 -127650    255300                         
## m2x    5 255248 255295 -127619    255238 62.779  1  2.313e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m3x <- lmer(finsat ~  religious_attendance + belief_god_centred  + age + sex + education_level + income_level + (1 | country), data = data, REML = FALSE)
anova(m2x,m3x, test = "Chisq")

## Data: data
## Models:
## m2x: finsat ~ religious_attendance + belief_god_centred + (1 | country)
## m3x: finsat ~ religious_attendance + belief_god_centred + age + sex + education_level + income_level + (1 | country)
##     npar    AIC    BIC  logLik -2*log(L) Chisq Df Pr(>Chisq)    
## m2x    5 255248 255295 -127619    255238                        
## m3x    9 244460 244545 -122221    244442 10796  4  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m4x <- lmer(finsat ~  religious_attendance + belief_god_centred  + age + sex +  education_level + income_level + hdi + mean_belief_god  + (1 | country), data = data, REML = FALSE)

anova(m3x,m4x, test = "Chisq")

## Data: data
## Models:
## m3x: finsat ~ religious_attendance + belief_god_centred + age + sex + education_level + income_level + (1 | country)
## m4x: finsat ~ religious_attendance + belief_god_centred + age + sex + education_level + income_level + hdi + mean_belief_god + (1 | country)
##     npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
## m3x    9 244460 244545 -122221    244442                     
## m4x   11 244461 244565 -122219    244439 2.8407  2     0.2416

screenreg(list(m1x, m2x, m3x, m4x), digits=3)

## 
## ============================================================================================
##                           Model 1          Model 2          Model 3          Model 4        
## --------------------------------------------------------------------------------------------
## (Intercept)                    -0.021           -0.020           -0.023           -0.017    
##                                (0.042)          (0.042)          (0.037)          (0.036)   
## religious_attendance            0.041 ***        0.030 ***        0.017 ***        0.017 ***
##                                (0.004)          (0.004)          (0.004)          (0.004)   
## belief_god_centred                               0.026 ***        0.037 ***        0.037 ***
##                                                 (0.003)          (0.003)          (0.003)   
## age                                                               0.038 ***        0.038 ***
##                                                                  (0.003)          (0.003)   
## sexMale                                                           0.005            0.005    
##                                                                  (0.006)          (0.006)   
## education_level                                                   0.025 ***        0.025 ***
##                                                                  (0.004)          (0.004)   
## income_level                                                      0.322 ***        0.322 ***
##                                                                  (0.003)          (0.003)   
## hdi                                                                                0.034    
##                                                                                   (0.053)   
## mean_belief_god                                                                   -0.035    
##                                                                                   (0.056)   
## --------------------------------------------------------------------------------------------
## AIC                        255308.326       255247.547       244459.649       244460.808    
## BIC                        255346.090       255294.753       244544.620       244564.661    
## Log Likelihood            -127650.163      -127618.773      -122220.824      -122219.404    
## Num. obs.                   93078            93078            93078            93078        
## Num. groups: country           63               63               63               63        
## Var: country (Intercept)        0.111            0.109            0.085            0.081    
## Var: Residual                   0.906            0.906            0.806            0.806    
## ============================================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05

model_a <- lmer(finsat ~  religious_attendance + belief_god_centred  + age + sex + education_level + income_level + hdi + mean_belief_god + belief_god_centred:hdi +  belief_god_centred:mean_belief_god + religious_attendance:mean_belief_god +  (1 + belief_god_centred + religious_attendance | country), data = data, REML = FALSE)

anova(m4x, model_a, test = "Chisq")

## Data: data
## Models:
## m4x: finsat ~ religious_attendance + belief_god_centred + age + sex + education_level + income_level + hdi + mean_belief_god + (1 | country)
## model_a: finsat ~ religious_attendance + belief_god_centred + age + sex + education_level + income_level + hdi + mean_belief_god + belief_god_centred:hdi + belief_god_centred:mean_belief_god + religious_attendance:mean_belief_god + (1 + belief_god_centred + religious_attendance | country)
##         npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
## m4x       11 244461 244565 -122219    244439                         
## model_a   19 244118 244297 -122040    244080 359.19  8  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_model(model_a)

	finsat
Predictors	Estimates	CI	p
(Intercept)	-0.01	-0.08 – 0.06	0.744
religious attendance	0.02	0.01 – 0.03	0.003
belief god centred	0.05	0.03 – 0.07	<0.001
age	0.04	0.03 – 0.05	<0.001
sex [Male]	0.00	-0.01 – 0.02	0.442
education level	0.02	0.02 – 0.03	<0.001
income level	0.32	0.32 – 0.33	<0.001
hdi	0.01	-0.09 – 0.11	0.815
mean belief god	-0.05	-0.16 – 0.05	0.333
belief god centred × hdi	-0.01	-0.04 – 0.02	0.385
belief god centred × mean belief god	0.02	-0.01 – 0.05	0.125
religious attendance × mean belief god	-0.01	-0.02 – 0.00	0.162
Random Effects
σ²	0.80
τ₀₀ _country	0.08
τ₁₁ _{country.belief_god_centred}	0.01
τ₁₁ _{country.religious_attendance}	0.00
ρ₀₁	0.55
	0.29
ICC	0.10
N _country	63
Observations	93078
Marginal R² / Conditional R²	0.117 / 0.203

model_b <- lmer(finsat ~  religious_attendance + belief_god_centred  + age + sex + education_level + income_level + mean_belief_god  +  belief_god_centred:mean_belief_god + religious_attendance:mean_belief_god +  (1 + belief_god_centred + religious_attendance | country), data = data, REML = FALSE)

tab_model(model_b)

	finsat
Predictors	Estimates	CI	p
(Intercept)	-0.01	-0.08 – 0.06	0.749
religious attendance	0.02	0.01 – 0.03	0.004
belief god centred	0.05	0.03 – 0.07	<0.001
age	0.04	0.03 – 0.05	<0.001
sex [Male]	0.00	-0.01 – 0.02	0.434
education level	0.02	0.02 – 0.03	<0.001
income level	0.32	0.32 – 0.33	<0.001
mean belief god	-0.06	-0.14 – 0.01	0.102
belief god centred × mean belief god	0.03	0.01 – 0.05	0.002
religious attendance × mean belief god	-0.01	-0.02 – 0.00	0.155
Random Effects
σ²	0.80
τ₀₀ _country	0.08
τ₁₁ _{country.belief_god_centred}	0.01
τ₁₁ _{country.religious_attendance}	0.00
ρ₀₁	0.54
	0.31
ICC	0.10
N _country	63
Observations	93078
Marginal R² / Conditional R²	0.116 / 0.202

model_c <- lmer(finsat ~  religious_attendance + belief_god_centred  + age + sex + education_level + income_level + hdi +  belief_god_centred:hdi  +  (1 + belief_god_centred + religious_attendance | country), data = data, REML = FALSE)

tab_model(model_a,model_b,model_c, show.ci = FALSE, title = "Table 2| Models Estimation", )

## Error in `tab_model()`:
## ! argument is missing, with no default

tab_model(m1x, m2x, m3x, m4x, model_a)

	finsat			finsat			finsat			finsat			finsat
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	-0.02	-0.10 – 0.06	0.618	-0.02	-0.10 – 0.06	0.624	-0.02	-0.10 – 0.05	0.542	-0.02	-0.09 – 0.05	0.641	-0.01	-0.08 – 0.06	0.744
religious attendance	0.04	0.03 – 0.05	<0.001	0.03	0.02 – 0.04	<0.001	0.02	0.01 – 0.02	<0.001	0.02	0.01 – 0.02	<0.001	0.02	0.01 – 0.03	0.003
belief god centred				0.03	0.02 – 0.03	<0.001	0.04	0.03 – 0.04	<0.001	0.04	0.03 – 0.04	<0.001	0.05	0.03 – 0.07	<0.001
age							0.04	0.03 – 0.04	<0.001	0.04	0.03 – 0.04	<0.001	0.04	0.03 – 0.05	<0.001
sex [Male]							0.00	-0.01 – 0.02	0.445	0.00	-0.01 – 0.02	0.443	0.00	-0.01 – 0.02	0.442
education level							0.03	0.02 – 0.03	<0.001	0.03	0.02 – 0.03	<0.001	0.02	0.02 – 0.03	<0.001
income level							0.32	0.32 – 0.33	<0.001	0.32	0.32 – 0.33	<0.001	0.32	0.32 – 0.33	<0.001
hdi										0.03	-0.07 – 0.14	0.526	0.01	-0.09 – 0.11	0.815
mean belief god										-0.03	-0.14 – 0.07	0.533	-0.05	-0.16 – 0.05	0.333
belief god centred × hdi													-0.01	-0.04 – 0.02	0.385
belief god centred × mean belief god													0.02	-0.01 – 0.05	0.125
religious attendance × mean belief god													-0.01	-0.02 – 0.00	0.162
Random Effects
σ²	0.91			0.91			0.81			0.81			0.80
τ₀₀	0.11 _country			0.11 _country			0.09 _country			0.08 _country			0.08 _country
τ₁₁													0.01 _{country.belief_god_centred}
													0.00 _{country.religious_attendance}
ρ₀₁													0.55
													0.29
ICC	0.11			0.11			0.10			0.09			0.10
N	63 _country			63 _country			63 _country			63 _country			63 _country
Observations	93078			93078			93078			93078			93078
Marginal R² / Conditional R²	0.002 / 0.110			0.002 / 0.109			0.109 / 0.194			0.117 / 0.198			0.117 / 0.203

We then add the control variables and analyse and interpret the final model.

Model Output

summary(model_a)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: finsat ~ religious_attendance + belief_god_centred + age + sex +  
##     education_level + income_level + hdi + mean_belief_god +  
##     belief_god_centred:hdi + belief_god_centred:mean_belief_god +  
##     religious_attendance:mean_belief_god + (1 + belief_god_centred +  
##     religious_attendance | country)
##    Data: data
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##  244117.6  244297.0 -122039.8  244079.6     93059 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6538 -0.6250  0.0294  0.6612  3.4575 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev. Corr     
##  country  (Intercept)          0.079844 0.28257           
##           belief_god_centred   0.005004 0.07074  0.55     
##           religious_attendance 0.001439 0.03793  0.29 0.08
##  Residual                      0.801978 0.89553           
## Number of obs: 93078, groups:  country, 63
## 
## Fixed effects:
##                                       Estimate Std. Error t value
## (Intercept)                          -0.011792   0.036113  -0.327
## religious_attendance                  0.018651   0.006344   2.940
## belief_god_centred                    0.051195   0.009876   5.184
## age                                   0.039157   0.003276  11.951
## sexMale                               0.004634   0.006029   0.769
## education_level                       0.023758   0.003538   6.716
## income_level                          0.322358   0.003180 101.368
## hdi                                   0.012105   0.051804   0.234
## mean_belief_god                      -0.053141   0.054870  -0.968
## belief_god_centred:hdi               -0.012739   0.014667  -0.869
## belief_god_centred:mean_belief_god    0.023210   0.015117   1.535
## religious_attendance:mean_belief_god -0.009441   0.006756  -1.397
## 
## Correlation of Fixed Effects:
##             (Intr) rlgs_t blf_g_ age    sexMal edctn_ incm_l hdi    mn_bl_
## rlgs_ttndnc  0.230                                                        
## blf_gd_cntr  0.492 -0.016                                                 
## age          0.003 -0.037 -0.019                                          
## sexMale     -0.077 -0.015  0.025 -0.052                                   
## educatn_lvl  0.002 -0.009  0.013  0.218 -0.046                            
## income_levl  0.000 -0.018  0.002  0.066 -0.023 -0.248                     
## hdi         -0.023  0.013 -0.006 -0.017  0.004 -0.017 -0.001              
## mean_blf_gd -0.089 -0.033 -0.036  0.007  0.000  0.004  0.002  0.720       
## blf_gd_cnt: -0.011 -0.016 -0.039 -0.006  0.010  0.010  0.002  0.503  0.364
## blf_gd_c:__ -0.040  0.012 -0.031 -0.001  0.009  0.006 -0.002  0.374  0.506
## rlgs_ttn:__ -0.058 -0.166  0.032 -0.004 -0.029  0.016  0.008  0.025  0.184
##             blf__: b__:__
## rlgs_ttndnc              
## blf_gd_cntr              
## age                      
## sexMale                  
## educatn_lvl              
## income_levl              
## hdi                      
## mean_blf_gd              
## blf_gd_cnt:              
## blf_gd_c:__  0.719       
## rlgs_ttn:__  0.003 -0.014

model_a <- lmer(finsat ~  religious_attendance + belief_god_centred  + age + sex + education_level + income_level + hdi + mean_belief_god + belief_god_centred:hdi +  belief_god_centred:mean_belief_god + religious_attendance:mean_belief_god +  (1 + belief_god_centred + religious_attendance | country), data = data)

tab_model(model_a)

	finsat
Predictors	Estimates	CI	p
(Intercept)	-0.01	-0.08 – 0.06	0.750
religious attendance	0.02	0.01 – 0.03	0.004
belief god centred	0.05	0.03 – 0.07	<0.001
age	0.04	0.03 – 0.05	<0.001
sex [Male]	0.00	-0.01 – 0.02	0.443
education level	0.02	0.02 – 0.03	<0.001
income level	0.32	0.32 – 0.33	<0.001
hdi	0.01	-0.09 – 0.12	0.818
mean belief god	-0.05	-0.16 – 0.06	0.346
belief god centred × hdi	-0.01	-0.04 – 0.02	0.393
belief god centred × mean belief god	0.02	-0.01 – 0.05	0.136
religious attendance × mean belief god	-0.01	-0.02 – 0.00	0.169
Random Effects
σ²	0.80
τ₀₀ _country	0.08
τ₁₁ _{country.belief_god_centred}	0.01
τ₁₁ _{country.religious_attendance}	0.00
ρ₀₁	0.55
	0.29
ICC	0.10
N _country	63
Observations	93078
Marginal R² / Conditional R²	0.116 / 0.206

effectsize::standardize_parameters(model_a, method = "pseudo")

## The following within-group terms have between-group variance:
##   religious_attendance, age, education_level, income_level
##   This can inflate standardized within-group parameters associated with
##   these terms.
##   See `help("demean", package = "datawizard")` for modeling between- and
##   within-subject effects.

## # Standardization method: pseudo
## 
## Parameter                              | Std. Coef. |        95% CI
## -------------------------------------------------------------------
## (Intercept)                            |       0.00 | [ 0.00, 0.00]
## religious attendance                   |       0.02 | [ 0.01, 0.03]
## belief god centred                     |       0.05 | [ 0.03, 0.07]
## age                                    |       0.04 | [ 0.03, 0.05]
## sex [Male]                             |   2.42e-03 | [ 0.00, 0.01]
## education level                        |       0.02 | [ 0.02, 0.03]
## income level                           |       0.34 | [ 0.33, 0.34]
## hdi                                    |       0.04 | [-0.28, 0.36]
## mean belief god                        |      -0.15 | [-0.47, 0.17]
## belief god centred × hdi               |      -0.01 | [-0.04, 0.02]
## belief god centred × mean belief god   |       0.02 | [-0.01, 0.06]
## religious attendance × mean belief god |  -9.29e-03 | [-0.02, 0.00]

It turns out that that the income level is the best indicator of individual financial satisfaction, followed by the country level of religiosity.

Second level analysis

dotplot(ranef(model_a, condVar=TRUE))

## $country

dotplot(coef(model_a))

## $country

As can be seen, our model shows that the random effect differs from one country to another. On average, Pakistanis have higher positive financial satisfaction, while Zimbabwe generally have negative financial satisfaction.

u <- ranef(model_a, condVar=TRUE)
u_df <- as.data.frame(u) # store as data frame for use with ggplot
u_df1 <- subset(u_df, term=="(Intercept)")
u_df2 <- subset(u_df, term=="religious_attendance")
u_df1$x <- u_df2$condval # u_df1 has intercept and slopes as separate columns
ggplot(u_df1, aes(x=condval, y=x)) + geom_point(shape=19, show.legend = FALSE) +
  geom_hline(yintercept=0) + geom_vline(xintercept=0) +
  labs(title="Figure B| Model A: Random Intercept vs Slope (Religious Attendance)", 
       caption="Dataset: WVS", x="Intercept", y="Slope")

u <- ranef(model_a, condVar=TRUE)
u_df <- as.data.frame(u) # store as data frame for use with ggplot
u_df1 <- subset(u_df, term=="(Intercept)")
u_df2 <- subset(u_df, term=="belief_god_centred")
u_df1$x <- u_df2$condval # u_df1 has intercept and slopes as separate columns
ggplot(u_df1, aes(x=condval, y=x)) + geom_point(shape=19, show.legend = FALSE) +
  geom_hline(yintercept=0) + geom_vline(xintercept=0) +
  labs(title="Figure C| Model A: Random Intercept vs Slope (Belief in God's Importance)", 
       caption="Dataset: WVS", x="Intercept", y="Slope")

mean_belief_vals <- quantile(data$mean_belief_god, probs = c(0.1, 0.5, 0.9), na.rm = TRUE)
belief_seq <- seq(min(data$belief_god_centred, na.rm = TRUE),
                  max(data$belief_god_centred, na.rm = TRUE), length.out = 50)

pred_grid <- expand.grid(
  belief_god_centred          = belief_seq,
  religious_attendance = mean(data$religious_attendance, na.rm = TRUE),
  age                  = mean(data$age, na.rm = TRUE),
  sex                  = names(which.max(table(data$sex))),
  income_level         = mean(data$income_level, na.rm = TRUE),
  education_level     = mean(data$education_level, na.rm = TRUE),
  mean_belief_god      = mean_belief_vals
) |>
  mutate(rel_label = factor(
    round(mean_belief_god, 2),
    labels = c("Low Religiosity (10th pct)", "Medium Religiosity (50th pct)", "High Religiosity (90th pct)")
  ))

# Predict using fixed effects only (re.form = NA ignores random effects)

pred_grid$finsat <- predict(model_b, newdata = pred_grid, re.form = NA)

# Plot
ggplot(pred_grid, aes(x = belief_god_centred, y = finsat, colour = rel_label, fill = rel_label)) +
  geom_line(linewidth = 1) +
  scale_colour_manual(
    values = c("Low Religiosity (10th pct)"    = "#E24B4A",
               "Medium Religiosity (50th pct)" = "#EF9F27",
               "High Religiosity (90th pct)"   = "#1D9E75"),
    name = "Religiosity level"
  ) +
  labs(
    title    = "Figure 1| Interaction: Belief in God × Country Religiosity",
    subtitle = "Predictions at mean values of all other covariates (fixed effects only)",
    x        = "Belief in God",
    y        = "Predicted Financial Satisfaction"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    legend.position  = "bottom",
    panel.grid.minor = element_blank(),
    plot.title       = element_text(face = "bold")
  )

interplot(model_c, "belief_god_centred", "hdi", hist=TRUE) + geom_hline(yintercept=0) +
    labs(title="Figure 2| Conditional Effect of Believing in God", caption="WVS round: 7",
         x="HDI (% of population)", y="Marginal Effect of believe in God")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the interplot package.
##   Please report the issue at <https://github.com/sammo3182/interplot/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the interplot package.
##   Please report the issue at <https://github.com/sammo3182/interplot/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Model Diagnostics

model_performance(model_a)

## Random effect variances not available. Returned R2 does not account for random effects.

## # Indices of model performance
## 
## AIC     |    AICc |     BIC | R2 (cond.) | R2 (marg.) |  RMSE | Sigma
## ---------------------------------------------------------------------
## 2.4e+05 | 2.4e+05 | 2.4e+05 |            |      0.128 | 0.895 | 0.896

The conditional R-squared value of 0.12 suggests that approximately 20% of the variance in the outcome variable can be explained by the fixed and random effects included in the model. This means that the predictors and factors considered in the model collectively account for 20% of the variability observed in the outcome variable. This is practically okay for a social science study based on Cohen (1988).

Linearity of residuals

#devtools::install_github("aloy/HLMdiag")
#install.packages("HLMdiag")

resid1_fm1 <- hlm_resid(model_b, level = 1, standardize = TRUE)
#resid1_fm1

ggplot(data = resid1_fm1, aes(x =religious_attendance, y = .std.ls.resid)) + 
  geom_point(alpha = 0.2) +
  geom_smooth(method = "loess", se = FALSE) + 
  labs(y = "LS level-1 residuals", 
       title = "LS residuals Religious Attendance")

## `geom_smooth()` using formula = 'y ~ x'

ggplot(data = resid1_fm1, aes(x =belief_god_centred, y = .std.ls.resid)) + 
  geom_point(alpha = 0.2) +
  geom_smooth(method = "loess", se = FALSE) + 
  labs(y = "LS level-1 residuals", 
       title = "belief in the importance of God")

## `geom_smooth()` using formula = 'y ~ x'

We can see that the smooth line runs horizontally through the residuals, suggesting that the residuals for the ‘belief in God’ variable are linearly related to the fitted values. Next we ran the Q-Q plot to analyse the linearity of residuals for the whole dataset.

qqnorm(residuals(model_b))
qqline(residuals(model_b))

Most of the residuals lie on the diagonal line, suggesting a linearity of our residuals. Therefore, our model confirm the assumption of the linearity of residuals.

Next we look at the second level residuals

resid1_fm2 <- hlm_resid(model_b, level = "country", standardize = TRUE)

## Warning: the 'nobars' function has moved to the reformulas package. Please update your imports, or ask an upstream package maintainter to do so.
## This warning is displayed once per session.

#head(resid1_fm2)

ggplot(data = resid1_fm2, aes(x = mean_belief_god, y = .std.ranef.intercept)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(y = "LS level-1 residuals",
       title = "Random against household size")

## `geom_smooth()` using formula = 'y ~ x'

As we can see, the smooth line divides the second-level residuals uniformly. This could confirm the linearity of the residuals. But we can spot an outlier in the top right-hand corner. Since this outlier is a country, removing it would harm our model. But let’s look at this outlier.

resid1_fm2[which(resid1_fm2$.std.ranef.intercept == max(resid1_fm2$.std.ranef.intercept)),]

## # A tibble: 1 × 8
##   country  mean_belief_god[,1] .std.ranef.intercept .std.ranef.belief_god_cent…¹
##   <chr>                  <dbl>                <dbl>                        <dbl>
## 1 Pakistan                1.02                 28.1                         4.81
## # ℹ abbreviated name: ¹.std.ranef.belief_god_centred
## # ℹ 4 more variables: .std.ranef.religious_attendance <dbl>,
## #   .std.ls.intercept <dbl>, .std.ls.belief_god_centred <dbl>,
## #   .std.ls.religious_attendance <dbl>

Pakistan is the outlier. On close examination of its distribution, we found nothing alarming to compel us to remove it from our dataset. We have therefore decided to keep it.

Homoscedasticity

Next we shall explore whether the variance of the model’s residuals is homogene across our predicted variable.

plot(fitted(model_b), residuals(model_b))

As our predictor is a quasi-interval variable, we obtain a strange-looking plot of the residuals against the fit. But an examination of the grouping of points across the bands allows us to confirm that the variance of the residuals is homogeneous. homoscedasticity assumption

Multicoleniarity

Next we look at the presence of Multicollinearity, among our IV’s.

vif(model_a)

##                 religious_attendance                   belief_god_centred 
##                             1.033569                             1.007660 
##                                  age                                  sex 
##                             1.071519                             1.006917 
##                      education_level                         income_level 
##                             1.136245                             1.084458 
##                                  hdi                      mean_belief_god 
##                             2.843851                             2.991550 
##               belief_god_centred:hdi   belief_god_centred:mean_belief_god 
##                             2.761873                             2.806694 
## religious_attendance:mean_belief_god 
##                             1.116258

Since none of the IVs has a VIF greater than 4, we should assert that there is no multicollineariy among our IV’s.

References

Berkessel, J. B., Gebauer, J. E., Joshanloo, M., Bleidorn, W., Rentfrow, P. J., Potter, J., & Gosling, S. D. (2021). National religiosity eases the psychological burden of poverty. Proceedings of the National Academy of Sciences, 118(39), e2103913118.
Cohen, J. (1988). The effect size. Statistical power analysis for the behavioral sciences, 77-83.
Cwynar, A., Potocki, T., Białowolski, P., & Węziak-Białowolska, D. (2024). Religious service attendance and consumer financial outcomes: Evidence from a longitudinal survey. Economics and Business Review, 10(4), 101-128. https://doi.org/10.18559/ebr.2024.4.1225
Friedman, H., & Lynch, J. (2023). Voluntary Simplicity as a Spiritual Remedy for Hypermaterialism and Overconsumption: Perspectives from Two of the Oldest Faiths. Journal of Intercultural Management and Ethics, 6(1), 55-70. https://doi.org/10.35478/jime.2023.1.07
Gebauer, J. E., Sedikides, C., Sch.nbrodt, F. D., Bleidorn, W., Rentfrow, P. J., Potter, J., & Gosling, S. D. (2017). The religiosity as social value hypothesis: A multi-method replication and extension across 65 countries and three levels of spatial aggregation. Journal of Personality and Social Psychology, 113(3), e18–e39. https://doi.org/10.1037/pspp0000104
Hong, S., Lee, J., Oh, F. D., & Shin, D. (2024). Religion and household saving behavior. Journal of Behavioral and Experimental Finance, 44, 100999.
Kalbasi, F., & Amani, M. (2022). Recognizing Extravagance Consumption and Avoiding Wastefulness Using Islamic Benchmarks in Nahj al-Balaghah. International Journal of Cultural & Religious Studies, 2(2). https://doi.org/10.32996/ijcrs.2022.2.2.1
Kortt, M., Dollery, B., & Grant, B. (2015). Religion and Life Satisfaction Down Under. Journal of Happiness Studies, 16, 277-293. https://doi.org/10.1007/S10902-014-9509-4.
Kose, T., & Cinar, K. (2024). A global assessment of the relationship between religiosity and financial satisfaction. The Social Science Journal, 61(2), 347–367. https://doi.org/10.1080/03623319.2020.1808769
Lim, C., & Putnam, R. (2010). Religion, Social Networks, and Life Satisfaction. American Sociological Review, 75, 914 - 933. https://doi.org/10.1177/0003122410386686.
Migheli, M. (2018). Religious Polarization, Religious Conflicts and Individual Financial Satisfaction: Evidence from India. Development Economics: Regional & Country Studies eJournal. https://doi.org/10.1111/rode.12567.
Mundlak, Y. (1978). On the pooling of time series and cross section data. Econometrica: journal of the Econometric Society, 69-85.

Sarofim, S., Minton, E., Hunting, A., Bartholomew, D., Zehra, S., Montford, W., Cabano, F., & Paul, P. (2020). Religion’s influence on the financial well-being of consumers: A conceptual framework and research agenda. Journal of Consumer Affairs, 54, 1028-1061. https://doi.org/10.1111/joca.12315.

Financial satisfaction and Religion

Paul Marie

2026-06-16

Financial satisfaction and Religion

Introduction

Exploratory Analysis

Variables selection

Variables Transformation

Evaluation and Description of Distribution

Missing Values and Imputation

Pattern of missingness

Imputation

Bivariates Tests

Financial satisfaction vs Numeric IVs

Financial satisfaction vs Categorical IVs

Multilevel regression

Forward Selection Models

Model Output

Model Diagnostics

Linearity of residuals

Homoscedasticity

Multicoleniarity

References