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.