library(readxl)
mydata <- read_excel("excelmagistrska.xlsx", sheet = "Sheet1", col_names = TRUE)
mydata <- mydata[-1, ] # drop the 2nd header row (full question text, not data)
head(mydata, 10)
## # A tibble: 10 × 41
## Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q11 Q12 Q13 Q14a Q14b Q14c Q14d Q15a Q15b
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 2 3 1 1 2 1 Petra pa… 3 -2 1 3 1 6 4 4 4 4 4 4
## 2 2 4 1 1 3 1 Parovel 2 -2 1 3 3 5 3 2 2 2 2 2
## 3 2 1 1 1 1 2 -2 -2 -2 -2 3 1 6 3 4 3 2 2 3
## 4 3 2 1 1 1 2 -2 -2 -2 -2 2 3 4 3 2 3 2 2 2
## 5 2 2 1 1 4 1 Se ne sp… 1 -2 2 3 4 6 4 4 4 2 3 4
## 6 2 1 1 1 2 1 Mojtribe 8 Craft 1 1 2 6 4 4 4 4 4 4
## 7 2 2 1 1 3 2 -2 -2 -2 -2 2 2 6 4 4 4 3 3 5
## 8 2 3 1 1 4 2 -2 -2 -2 -2 1 4 6 2 3 3 3 1 2
## 9 4 2 1 1 3 2 -2 -2 -2 -2 2 2 5 3 4 4 3 4 4
## 10 3 3 1 1 3 1 Aeon 8 Aeon 1 3 2 3 3 3 3 3 3 3
## # ℹ 22 more variables: Q15c <chr>, Q15d <chr>, Q16a <chr>, Q16b <chr>, Q16c <chr>, Q16d <chr>, Q17a <chr>, Q17b <chr>,
## # Q17c <chr>, Q17d <chr>, Q18a <chr>, Q18b <chr>, Q18c <chr>, Q18d <chr>, Q19a <chr>, Q19b <chr>, Q19c <chr>,
## # Q19d <chr>, Q20 <chr>, Q21 <chr>, Q22 <chr>, Q23 <chr>
#Identify values -1,-2 and replace them with NA. -3 values were deleted in excel
# -1 means: The person has not answered the specific question
# -2 means: Person did not respond because it did not satisfy if sentence
mydata[mydata == -1] <- NA
mydata[mydata == -2] <- NA
num_cols <- setdiff(names(mydata), c("Q7", "Q9"))
mydata[num_cols] <- lapply(mydata[num_cols], function(x) suppressWarnings(as.numeric(x)))
head(mydata, 10)
## # A tibble: 10 × 41
## Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q11 Q12 Q13 Q14a Q14b Q14c Q14d Q15a Q15b
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 3 1 1 2 1 Petra pa… 3 <NA> 1 3 1 6 4 4 4 4 4 4
## 2 2 4 1 1 3 1 Parovel 2 <NA> 1 3 3 5 3 2 2 2 2 2
## 3 2 1 1 1 1 2 <NA> NA <NA> NA 3 1 6 3 4 3 2 2 3
## 4 3 2 1 1 1 2 <NA> NA <NA> NA 2 3 4 3 2 3 2 2 2
## 5 2 2 1 1 4 1 Se ne sp… 1 <NA> 2 3 4 6 4 4 4 2 3 4
## 6 2 1 1 1 2 1 Mojtribe 8 Craft 1 1 2 6 4 4 4 4 4 4
## 7 2 2 1 1 3 2 <NA> NA <NA> NA 2 2 6 4 4 4 3 3 5
## 8 2 3 1 1 4 2 <NA> NA <NA> NA 1 4 6 2 3 3 3 1 2
## 9 4 2 1 1 3 2 <NA> NA <NA> NA 2 2 5 3 4 4 3 4 4
## 10 3 3 1 1 3 1 Aeon 8 Aeon 1 3 2 3 3 3 3 3 3 3
## # ℹ 22 more variables: Q15c <dbl>, Q15d <dbl>, Q16a <dbl>, Q16b <dbl>, Q16c <dbl>, Q16d <dbl>, Q17a <dbl>, Q17b <dbl>,
## # Q17c <dbl>, Q17d <dbl>, Q18a <dbl>, Q18b <dbl>, Q18c <dbl>, Q18d <dbl>, Q19a <dbl>, Q19b <dbl>, Q19c <dbl>,
## # Q19d <dbl>, Q20 <dbl>, Q21 <dbl>, Q22 <dbl>, Q23 <dbl>
Description:
mydata$GenderF <- factor(mydata$Q20,
levels = c(1, 2),
labels = c("Male", "Female"))
mydata$EducationF <- factor(mydata$Q21,
levels = c(1, 2, 3, 4, 5, 6),
labels = c("Primary school", "Lower or vocational secondary school", "Upper secondary", "Higher vocational education", "University degree", "Postgraduate"))
mydata$AreaF <- factor(mydata$Q22,
levels = c(1, 2, 3),
labels = c("Urban", "Suburban", "Rural"))
mydata$Credibility <- rowMeans(mydata[, c("Q14a", "Q14b", "Q14c", "Q14d")], na.rm = FALSE)
mydata$Authenticity <- rowMeans(mydata[, c("Q15a", "Q15b", "Q15c", "Q15d")], na.rm = FALSE)
mydata$Engagement <- rowMeans(mydata[, c("Q16a", "Q16b", "Q16c", "Q16d")], na.rm = FALSE)
mydata$Parasocial <- rowMeans(mydata[, c("Q17a", "Q17b", "Q17c", "Q17d")], na.rm = FALSE)
mydata$Homophily <- rowMeans(mydata[, c("Q18a", "Q18b", "Q18c", "Q18d")], na.rm = FALSE)
mydata$BrandLoyalty <- rowMeans(mydata[, c("Q19a", "Q19b", "Q19c", "Q19d")], na.rm = FALSE)
model_vars <- c("BrandLoyalty", "Credibility", "Authenticity", "Engagement",
"Parasocial", "Homophily")
mydata <- mydata[complete.cases(mydata[, model_vars]), ]
row.names(mydata) <- NULL
library(pastecs)
#descriptive stat. on the constructs
round(stat.desc(mydata[, c("BrandLoyalty", "Credibility", "Authenticity",
"Engagement", "Parasocial", "Homophily")]), 2)
## BrandLoyalty Credibility Authenticity Engagement Parasocial Homophily
## nbr.val 73.00 73.00 73.00 73.00 73.00 73.00
## nbr.null 0.00 0.00 0.00 0.00 0.00 0.00
## nbr.na 0.00 0.00 0.00 0.00 0.00 0.00
## min 2.00 2.00 1.75 1.00 1.00 1.00
## max 5.00 5.00 5.00 4.75 5.00 5.00
## range 3.00 3.00 3.25 3.75 4.00 4.00
## sum 239.75 273.50 279.50 203.00 238.00 243.00
## median 3.25 3.75 4.00 2.50 3.50 3.25
## mean 3.28 3.75 3.83 2.78 3.26 3.33
## SE.mean 0.10 0.09 0.09 0.11 0.11 0.10
## CI.mean.0.95 0.19 0.18 0.19 0.22 0.23 0.19
## var 0.69 0.57 0.64 0.93 0.96 0.67
## std.dev 0.83 0.75 0.80 0.96 0.98 0.82
## coef.var 0.25 0.20 0.21 0.35 0.30 0.25
scatterplotMatrix(mydata[, c("BrandLoyalty", "Credibility", "Authenticity",
"Engagement", "Parasocial", "Homophily")],
smooth = FALSE)
library(Hmisc)
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
rcorr(as.matrix(mydata[, c("BrandLoyalty", "Credibility", "Authenticity",
"Engagement", "Parasocial", "Homophily")]))
## BrandLoyalty Credibility Authenticity Engagement Parasocial Homophily
## BrandLoyalty 1.00 0.58 0.47 0.67 0.66 0.49
## Credibility 0.58 1.00 0.82 0.55 0.66 0.53
## Authenticity 0.47 0.82 1.00 0.54 0.67 0.45
## Engagement 0.67 0.55 0.54 1.00 0.75 0.51
## Parasocial 0.66 0.66 0.67 0.75 1.00 0.71
## Homophily 0.49 0.53 0.45 0.51 0.71 1.00
##
## n= 73
##
##
## P
## BrandLoyalty Credibility Authenticity Engagement Parasocial Homophily
## BrandLoyalty 0 0 0 0 0
## Credibility 0 0 0 0 0
## Authenticity 0 0 0 0 0
## Engagement 0 0 0 0 0
## Parasocial 0 0 0 0 0
## Homophily 0 0 0 0 0
A Pearson correlation coefficient was computed to assess the relationship between Authenticity and Credibility. There was a strong positive correlation between the two variables.
DIAGNOSTIC MODEL - just to check assumptions before we interpret anything
fit1 <- lm(BrandLoyalty ~ Credibility + Authenticity + Engagement +
Parasocial + Homophily,
data = mydata)
summary(fit1)
##
## Call:
## lm(formula = BrandLoyalty ~ Credibility + Authenticity + Engagement +
## Parasocial + Homophily, data = mydata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.27056 -0.44078 0.01961 0.42281 1.17126
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.968487 0.386511 2.506 0.01466 *
## Credibility 0.416796 0.167615 2.487 0.01539 *
## Authenticity -0.249102 0.157295 -1.584 0.11798
## Engagement 0.333882 0.106669 3.130 0.00259 **
## Parasocial 0.231911 0.139628 1.661 0.10140
## Homophily 0.007028 0.120363 0.058 0.95361
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5752 on 67 degrees of freedom
## Multiple R-squared: 0.5514, Adjusted R-squared: 0.5179
## F-statistic: 16.47 on 5 and 67 DF, p-value: 1.424e-10
#Multicollinearity (want close to 1, definitely not above 5)
vif(fit1)
## Credibility Authenticity Engagement Parasocial Homophily
## 3.454675 3.457208 2.290866 4.090931 2.106444
mean(vif(fit1))
## [1] 3.080025
#standardized residuals - outliers, +Cook’s distance - units with large influence
mydata$StdResid <- round(rstandard(fit1), 3)
mydata$CooksD <- round(cooks.distance(fit1), 3)
hist(mydata$StdResid,
xlab = "Standardized residuals",
ylab = "Frequency",
main = "Histogram of standardized residuals")
shapiro.test(mydata$StdResid)
##
## Shapiro-Wilk normality test
##
## data: mydata$StdResid
## W = 0.98627, p-value = 0.6153
p < .05 → you reject H0 → IS NORMALLY DISTRIBUTED p > .05 → you fail to reject H0 → IS NOT NORMALLY DISTRIBUTED
hist(mydata$CooksD,
xlab = "Cooks distance",
ylab = "Frequency",
main = "Histogram of Cooks distances")
We dont have any distances above 1.
head(mydata[order(mydata$StdResid), c("StdResid", "CooksD")], 5)
## # A tibble: 5 × 2
## StdResid CooksD
## <dbl> <dbl>
## 1 -2.25 0.03
## 2 -1.75 0.065
## 3 -1.74 0.113
## 4 -1.70 0.025
## 5 -1.65 0.021
head(mydata[order(-mydata$CooksD), c("StdResid", "CooksD")], 5)
## # A tibble: 5 × 2
## StdResid CooksD
## <dbl> <dbl>
## 1 1.67 0.123
## 2 2.18 0.121
## 3 2.18 0.114
## 4 -1.74 0.113
## 5 2.02 0.1
mydata$StdFittedValues <- scale(fit1$fitted.values)
scatterplot(y = mydata$StdResid, x = mydata$StdFittedValues,
ylab = "Standardized residuals",
xlab = "Standardized fitted values",
boxplots = FALSE,
regLine = FALSE,
smooth = FALSE)
Linearity OK points look random around 0 across the whole X range.
Homoscedasticity OK since vertical spread of points looks roughly
constant left to right (no funnel/fan shape).
#formal test, H0: variance is constant (homoscedastic)
library(olsrr)
##
## Attaching package: 'olsrr'
## The following object is masked from 'package:datasets':
##
## rivers
ols_test_breusch_pagan(fit1)
##
## Breusch Pagan Test for Heteroskedasticity
## -----------------------------------------
## Ho: the variance is constant
## Ha: the variance is not constant
##
## Data
## ----------------------------------------
## Response : BrandLoyalty
## Variables: fitted values of BrandLoyalty
##
## Test Summary
## ----------------------------
## DF = 1
## Chi2 = 0.1017924
## Prob > Chi2 = 0.7496892
The Breusch–Pagan test was conducted to assess the assumption of homoscedasticity. The test was not statistically significant, χ²(1) = 0.10, p = .750, indicating no evidence of heteroskedasticity. Therefore, the assumption of constant error variance was considered satisfied.
summary(fit1)
##
## Call:
## lm(formula = BrandLoyalty ~ Credibility + Authenticity + Engagement +
## Parasocial + Homophily, data = mydata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.27056 -0.44078 0.01961 0.42281 1.17126
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.968487 0.386511 2.506 0.01466 *
## Credibility 0.416796 0.167615 2.487 0.01539 *
## Authenticity -0.249102 0.157295 -1.584 0.11798
## Engagement 0.333882 0.106669 3.130 0.00259 **
## Parasocial 0.231911 0.139628 1.661 0.10140
## Homophily 0.007028 0.120363 0.058 0.95361
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5752 on 67 degrees of freedom
## Multiple R-squared: 0.5514, Adjusted R-squared: 0.5179
## F-statistic: 16.47 on 5 and 67 DF, p-value: 1.424e-10
CREDIBILITY: H0: β1 = 0 H1: β1 ≠ 0 p<0.05. We can reject null hypothesis. With every additional increase in one unit of credibility, an average brand loyalty will increase by 0.41, assuming that everything else (authenticity, engagement, parasocial relationship and homophily) remains unchanged.
AUTHENTICITY: H0: β2 = 0 H1: β2 ≠ 0 p>0.05. We cannot reject null hypothesis. We did not find significant effect of Authenticity on Brand Loyalty among Slovenian Gen Z population, assuming that everything else (credibility, engagement, parasocial relationship and homophily) remains unchanged.
ENGAGEMENT: H0: β3 = 0 H1: β3 ≠ 0 p<0.05. We can reject null hypothesis. With every additional increase in one unit of engagement, an average brand loyalty will increase by 0.33, assuming that everything else (authenticity, credibility, parasocial relationship and homophily) remains unchanged.
PARASOCIAL RELATIONSHIP: H0: β4 = 0 H1: β4 ≠ 0 p>0.05. We cannot reject null hypothesis. We did not find significant effect of Parasocial relationships on Brand Loyalty among Slovenian Gen Z population, assuming that everything else (credibility, authenticity, engagement and homophily) remains unchanged.
HOMOPHILY: H0: β5 = 0 H1: β5 ≠ 0 p>0.05. We cannot reject null hypothesis. We did not find significant effect of homophily on Brand Loyalty among Slovenian Gen Z population, assuming that everything else (credibility, authenticity, engagement and parasocial relationship) remains unchanged.
R squared: 0.55; 55% of variability of brand loyalty is explained by linear effect of authenticity, engagement, parasocial relationship and homophily.
ANOVA: H0: RO2=0 H1: RO2>0 We can reject Ho at p<0,001. We found the linear variability between dependent and at least one explanatory variables.
sqrt(summary(fit1)$r.squared)
## [1] 0.7425414
Multiple coefficient of correlation: 0,742 - linear relationship between brand loyalty and all 5 explanatory variables is strong.
#standardized coefficients - which variabvle contributes most? (H2-H5 relative importance)
library(lm.beta)
lm.beta(fit1)
##
## Call:
## lm(formula = BrandLoyalty ~ Credibility + Authenticity + Engagement +
## Parasocial + Homophily, data = mydata)
##
## Standardized Coefficients::
## (Intercept) Credibility Authenticity Engagement Parasocial Homophily
## NA 0.378201014 -0.240954005 0.387670889 0.274895442 0.006934678
Highest are credibility and engagement, which is expected since they were significant.