Corporate responsibility Versus Easiness of E-Commerce Payments Methods: An Experimental Study

1. Introduction

This case study examines the impact of introducing a new payment method, “ControvPay,” on customer perceptions. While the addition was anticipated to enhance convenience, it also risked provoking negative reactions due to ethical controversies associated with the provider. To evaluate this trade-off, we conducted an experimental study measuring key metrics like adoption, perceived ease of use, perceived social responsability, purchase intention, and loyalty.

The analytical approach encompasses:

• Data simulation and cleaning.
• Reliability analysis.
• Exploratory analysis with descriptive statistics and visualizations.
• Balance checks using inferential tests (chi-square, t-test, ANOVA, correlations).
• Comparative analyses with factorial ANOVA and ANCOVA models.
• Sensitivity analyses with Bayesian modeling.
• Mediation Analyses.
• Simulation-based predictions.

Although inspired by a real-world business challenge, the results presented here are simulated to ensure confidentiality.

2. Setup

We begin by loading the R packages required for data manipulation, visualization, and statistical modeling.

# Data handling
library(dplyr)       # Data manipulation

## 
## Adjuntando el paquete: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# Visualization
library(ggplot2)     # General plotting
library(kableExtra)  # Table formatting

## 
## Adjuntando el paquete: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

library(DiagrammeR)  # Path Diagrams 

# Statistical analysis
library(psych)       # Descriptive statistics

## 
## Adjuntando el paquete: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(MASS)        # Simulation of multivariate distributions

## 
## Adjuntando el paquete: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(car)         # Scatterplot matrices

## Cargando paquete requerido: carData

## 
## Adjuntando el paquete: 'car'

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

library(corrplot)    # Correlation plots

## corrplot 0.95 loaded

library(effectsize)  # Effect size calculations

## 
## Adjuntando el paquete: 'effectsize'

## The following object is masked from 'package:psych':
## 
##     phi

library(BayesFactor) # Bayesian inference

## Cargando paquete requerido: coda

## Cargando paquete requerido: Matrix

## ************
## Welcome to BayesFactor 0.9.12-4.7. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
## 
## Type BFManual() to open the manual.
## ************

# Mediation analysis
# For mediation analyses we used PROCESS for R, available at: https://www.processmacro.org/index.html

3. Data Simulation

3.1. Methodological Brief

Data for this study were simulated to replicate a survey experiment where participants were randomly assigned to one of two conditions. In the experimental condition, participants viewed an e-commerce platform interface that included the controversial payment method “ControvPay” alongside standard options (e.g., credit card, bank transfer…). The control condition presented the identical platform but without ControvPay. After exposure to their assigned condition, participants completed questionnaires assessing their perceptions and demographic characteristics.

3.2. Experimental Assignment

We initialized the data simulation by randomly assigning participants to two experimental groups:

• Control: Viewed the payment page with standard payment options only.
• Experimental: Viewed the same payment page with the addition of ControvPay.

# Set seed for reproducibility
set.seed(123)

# Define dataset size
data.n <- 800

# Initialize dataframe with IDs
pay.df <- data.frame(id = factor(c(1:data.n)))

# Create the factor about the experimental assignment
pay.df$condition <- factor(sample(x = c("control", "experimental"), size = data.n, replace = TRUE))

# Check the data distributions
summary(pay.df)

##        id             condition  
##  1      :  1   control     :402  
##  2      :  1   experimental:398  
##  3      :  1                     
##  4      :  1                     
##  5      :  1                     
##  6      :  1                     
##  (Other):794

3.3. Participants’ Information

In addition to basic demographic variables of age and gender, we simulated purchasing frequency at the retailer during the previous year. This question was embedded among similar questions about other retailers to maintain anonymity and provide contextual framing for the survey.

# Simulate age with log-normal distribution, capped between 15 and 85
pay.df$age <- round(rlnorm(n = data.n, meanlog = log(40), sdlog = log(1.4)))
pay.df$age[pay.df$age > 85] <- 85
pay.df$age[pay.df$age < 15] <- 15

# Simulate gender
pay.df$gender <- factor(sample(x = c("man", "woman"),
                              size = data.n,
                              replace = TRUE,
                              prob = c(0.5, 0.5)))

# Simulate purchasing frequency at the retailer during the last year
pay.df$freq.customer <- factor(sample(x = c("none", "1 time", "2-4 times", "5-9 times", "10 times or more"),
                                      size = data.n,
                                      replace = TRUE,
                                      prob = c(0.5, 0.2, 0.15, 0.10, 0.05)),
                               ordered = TRUE, 
                               levels = c("none", "1 time", "2-4 times", "5-9 times", "10 times or more"))

# Create binary customer indicator based on purchasing frequency
pay.df$customer <- rep(NA_character_, data.n)
pay.df$customer[pay.df$freq.customer == "none"] <- "no"
pay.df$customer[pay.df$freq.customer %in% c("1 time", "2-4 times", "5-9 times", "10 times or more")] <- "yes"
pay.df$customer <- factor(pay.df$customer)

# Check data distributions
summary(pay.df)

##        id             condition        age          gender   
##  1      :  1   control     :402   Min.   :16.00   man  :399  
##  2      :  1   experimental:398   1st Qu.:32.00   woman:401  
##  3      :  1                      Median :40.50              
##  4      :  1                      Mean   :42.53              
##  5      :  1                      3rd Qu.:50.00              
##  6      :  1                      Max.   :85.00              
##  (Other):794                                                 
##           freq.customer customer 
##  none            :396   no :396  
##  1 time          :162   yes:404  
##  2-4 times       :110            
##  5-9 times       : 89            
##  10 times or more: 43            
##                                  
## 

3.4. Questionnaire Responses

For each dependent variable, we present the questionnaire items and describe their data simulation approach.

For multi-item questionnaires, we simulated Likert-type responses using a custom function that employs multivariate normal distributions from the MASS package, as specified below.

# Function for simulating Likert-type items in multi-item questionnaires

# Function name: f.sim.likert
#
# Input parameters:
#   n: number of cases to simulate
#   n.items: number of items to simulate  
#   mean: target mean for the items
#   sd: target standard deviation for the items
#   rounding: if TRUE, rounds scores and constrains within scale limits (default: TRUE)
#   scale.max: maximum score of the Likert scale
#   scale.min: minimum score of the Likert scale
#   corr: inter-item correlation (default: 0.6)
#
# Output: simulated Likert-type item scores using multivariate normal distribution

f.sim.likert <- function(n, n.items, mean, sd, rounding = TRUE, scale.max, scale.min, corr = 0.6) {

  # Ensure MASS package is available
  if (!requireNamespace("MASS", quietly = TRUE)) {
    install.packages("MASS")
  }
  if (!"package:MASS" %in% search()) {
    library(MASS, quietly = TRUE)
  }
  
  # Specify correlation matrix
  cor.matrix <- matrix(corr, nrow = n.items, ncol = n.items)
  diag(cor.matrix) <- 1
  
  # Calculate covariance matrix
  cov.matrix <- cor.matrix * (sd^2)
  diag(cov.matrix) <- sd^2
  
  # Simulate correlated scores
  scores <- mvrnorm(n, mu = rep(mean, n.items), Sigma = cov.matrix)
  
  # Apply rounding and scale constraints if requested
  if(rounding == TRUE) {
    scores <- round(pmin(pmax(scores, scale.min), scale.max))
  }
  return(scores)
}

We now apply this function to simulate responses for the multi-item scales.

We now apply this function to simulate responses for the multi-item scales.

3.4.1. Perceived Corporate Social Responsibility

The questionnaire items measuring perceptions of the website’s social responsibility are presented below.

# Dataframe to display items
measure.sr <- data.frame(
  Item = c("CSR1", "CSR2", "CSR3", "CSR4"),
  Text = c(
    "I believe that this retailer ensures that the respect of ethical principles has priority over economic performance.",
    "I believe that this retailer permits ethical concerns to negatively affect economic performance.",
    "I believe that this retailer is committed to well-defined ethics principles.",
    "I believe that this retailer avoids compromising ethical standards in order to achieve corporate goals."
  )
) 

# Create formatted table with measure description
kable(measure.sr, "html", booktabs = TRUE,
      caption = "Table 1. Perception of Social Responsibility") %>%
  add_footnote(
    label = "Note: Items were adapted from Maignan (2001). Participants responded using a 7-point Likert scale from 1 (strongly disagree) to 7 (completely agree).",
    notation = "none"
  )

Table 1. Perception of Social Responsibility

Item	Text
CSR1	I believe that this retailer ensures that the respect of ethical principles has priority over economic performance.
CSR2	I believe that this retailer permits ethical concerns to negatively affect economic performance.
CSR3	I believe that this retailer is committed to well-defined ethics principles.
CSR4	I believe that this retailer avoids compromising ethical standards in order to achieve corporate goals.
Note: Items were adapted from Maignan (2001). Participants responded using a 7-point Likert scale from 1 (strongly disagree) to 7 (completely agree).

We simulated responses to the social responsibility items with slightly lower scores in the experimental condition, but only among the retailer’s existing customers.

# Simulate preliminary distribution using function from section 3.4
csr.scores <- f.sim.likert(n = data.n, 
                           n.items = 4,
                           mean = 5.5, 
                           sd = 1.5,
                           rounding = FALSE,
                           scale.max = 7,
                           scale.min = 1, 
                           corr = 0.6)

# Apply experimental effect: lower scores for customers in experimental condition
effect.size <- 0.35 
csr.scores[pay.df$condition == "experimental" & pay.df$customer == "yes"] <- 
  csr.scores[pay.df$condition == "experimental" & pay.df$customer == "yes"] - effect.size

# Round scores and constrain within scale limits
csr.scores <- round(pmin(pmax(csr.scores, 1), 7))

# Assign scores to dataframe
pay.df$csr1 <- csr.scores[, 1]
pay.df$csr2 <- csr.scores[, 2]
pay.df$csr3 <- csr.scores[, 3]
pay.df$csr4 <- csr.scores[, 4]

# Check distributions and relationships
par(mfrow = c(2, 2))
barplot(table(pay.df$csr1), main = "CSR1 Distribution")
barplot(table(pay.df$csr2), main = "CSR2 Distribution") 
barplot(table(pay.df$csr3), main = "CSR3 Distribution")
barplot(table(pay.df$csr4), main = "CSR4 Distribution")

# Check inter-item correlations
corr.test(pay.df[, c("csr1", "csr2", "csr3", "csr4")])

## Call:corr.test(x = pay.df[, c("csr1", "csr2", "csr3", "csr4")])
## Correlation matrix 
##      csr1 csr2 csr3 csr4
## csr1 1.00 0.52 0.55 0.54
## csr2 0.52 1.00 0.54 0.55
## csr3 0.55 0.54 1.00 0.53
## csr4 0.54 0.55 0.53 1.00
## Sample Size 
## [1] 800
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##      csr1 csr2 csr3 csr4
## csr1    0    0    0    0
## csr2    0    0    0    0
## csr3    0    0    0    0
## csr4    0    0    0    0
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

# Examine scores across experimental conditions
aggregate(pay.df$csr1, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

aggregate(pay.df$csr2, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

aggregate(pay.df$csr3, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

aggregate(pay.df$csr4, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

The simulation performed as expected, showing slightly lower social responsibility perceptions among existing customers exposed to ControvPay.

3.4.2. Perceived Ease of Use

The questionnaire items measuring expected ease of the payment process are presented below.

# Dataframe to display items
measure.eou <- data.frame(
  Item = c("EOU1", "EOU2", "EOU3", "EOU4"),
  Text = c(
    "My interaction with the payment service will be clear and understandable.",
    "Interacting with the payment service will not require a lot of my mental effort.",
    "I expect to find the payment service easy to use.",
    "I expect to find it easy to get the payment service to do what I want it to do."
  )
) 

# Create formatted table
kable(measure.eou, format = "html", 
      caption = "Table 2. Perceived Ease of Use") %>%
  add_footnote(
    label = "Note: Items were adapted from the Ease of Use scale by Venkatesh & Davis (2000). Participants responded using a 7-point Likert scale from 1 (strongly disagree) to 7 (completely agree).",
    notation = "none"
  )

Table 2. Perceived Ease of Use

Item	Text
EOU1	My interaction with the payment service will be clear and understandable.
EOU2	Interacting with the payment service will not require a lot of my mental effort.
EOU3	I expect to find the payment service easy to use.
EOU4	I expect to find it easy to get the payment service to do what I want it to do.
Note: Items were adapted from the Ease of Use scale by Venkatesh & Davis (2000). Participants responded using a 7-point Likert scale from 1 (strongly disagree) to 7 (completely agree).

We simulated ease of use responses as independent of experimental condition but inversely correlated with participant age.

# Simulate preliminary distribution using function from section 3.4
eou.scores <- f.sim.likert(n = data.n, 
                           n.items = 4,
                           mean = 6, 
                           sd = 1.6,
                           rounding = FALSE,
                           corr = 0.6)

# Apply age effect: lower scores for older participants
age.effect.size <- 0.04 

for(i in 1:ncol(eou.scores)) {
  eou.scores[, i] <- eou.scores[, i] - pay.df$age * age.effect.size
}

# Ensure scores remain within scale limits
eou.scores <- round(pmin(pmax(eou.scores, 1), 7))

# Assign scores to dataframe
pay.df$eou1 <- eou.scores[, 1]
pay.df$eou2 <- eou.scores[, 2]
pay.df$eou3 <- eou.scores[, 3]
pay.df$eou4 <- eou.scores[, 4]

# Check distributions and relationships
par(mfrow = c(2, 2))
barplot(table(pay.df$eou1), main = "EOU1 Distribution")
barplot(table(pay.df$eou2), main = "EOU2 Distribution")
barplot(table(pay.df$eou3), main = "EOU3 Distribution") 
barplot(table(pay.df$eou4), main = "EOU4 Distribution")

# Check inter-item correlations
corr.test(pay.df[, c("eou1", "eou2", "eou3", "eou4")])

## Call:corr.test(x = pay.df[, c("eou1", "eou2", "eou3", "eou4")])
## Correlation matrix 
##      eou1 eou2 eou3 eou4
## eou1 1.00 0.61 0.63 0.62
## eou2 0.61 1.00 0.63 0.61
## eou3 0.63 0.63 1.00 0.61
## eou4 0.62 0.61 0.61 1.00
## Sample Size 
## [1] 800
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##      eou1 eou2 eou3 eou4
## eou1    0    0    0    0
## eou2    0    0    0    0
## eou3    0    0    0    0
## eou4    0    0    0    0
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

# Examine scores across experimental conditions
aggregate(pay.df$eou1, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

aggregate(pay.df$eou2, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

aggregate(pay.df$eou3, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

aggregate(pay.df$eou4, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

3.4.3. Purchase Intention

The questionnaire items measuring purchase intention after viewing the payment options are presented below.

# Dataframe to display items
measure.pi <- data.frame(
  Item = c("PI1", "PI2"),
  Text = c(
    "I would intend to use the payment options.",
    "I predict that I would use the payment options."
  )
) 

# Create formatted table
kable(measure.pi, format = "html", 
      caption = "Table 3. Purchase Intention") %>%
  add_footnote(
    label = "Note: Items were adapted from the Intention to Use scale by Venkatesh & Davis (2000). Participants responded using a 7-point Likert scale from 1 (strongly disagree) to 7 (completely agree).",
    notation = "none"
  )

Table 3. Purchase Intention

Item	Text
PI1	I would intend to use the payment options.
PI2	I predict that I would use the payment options.
Note: Items were adapted from the Intention to Use scale by Venkatesh & Davis (2000). Participants responded using a 7-point Likert scale from 1 (strongly disagree) to 7 (completely agree).

We simulated purchase intention responses as moderately influenced by perceived ease of use and slightly influenced by social responsibility perceptions.

# Simulate preliminary distribution using function from section 3.4
pi.scores <- f.sim.likert(n = data.n, 
                          n.items = 2,
                          mean = 5, 
                          sd = 1.2,
                          rounding = FALSE,
                          scale.max = 7,
                          scale.min = 1, 
                          corr = 0.6)

# Apply relationships with ease of use and social responsibility perceptions

# Calculate standardized scale means
csr.mean <- scale(rowMeans(pay.df[, c("csr1", "csr2", "csr3", "csr4")]))
eou.mean <- scale(rowMeans(pay.df[, c("eou1", "eou2", "eou3", "eou4")]))

# Define effect sizes
csr.effect.size <- 0.1
eou.effect.size <- 0.3 

# Apply effects to each item
for(i in 1:ncol(pi.scores)) {
  pi.scores[, i] <- pi.scores[, i] + csr.mean * csr.effect.size + eou.mean * eou.effect.size
}

# Ensure scores remain within scale limits
pi.scores <- round(pmin(pmax(pi.scores, 1), 7))

# Assign scores to dataframe
pay.df$pi1 <- pi.scores[, 1]
pay.df$pi2 <- pi.scores[, 2]

# Check distributions and relationships
par(mfrow = c(1, 2))
barplot(table(pay.df$pi1), main = "PI1 Distribution")
barplot(table(pay.df$pi2), main = "PI2 Distribution")

# Check inter-item correlations
corr.test(pay.df[, c("pi1", "pi2")])

## Call:corr.test(x = pay.df[, c("pi1", "pi2")])
## Correlation matrix 
##     pi1 pi2
## pi1 1.0 0.6
## pi2 0.6 1.0
## Sample Size 
## [1] 800
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##     pi1 pi2
## pi1   0   0
## pi2   0   0
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

# Examine scores across experimental conditions
aggregate(pay.df$pi1, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

aggregate(pay.df$pi2, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

3.4.4. Customer Loyalty

The questionnaire items measuring customer loyalty after viewing the payment options are presented below.

# Dataframe to display items
measure.cl <- data.frame(
  Item = c("CL1", "CL2"),
  Text = c(
    "I will likely visit the website in the future.",
    "How likely are you to recommend the website to a friend or colleague?"
  )
)

# Create formatted table
kable(measure.cl, format = "html", 
      caption = "Table 4. Customer Loyalty") %>%
  add_footnote(
    label = "Note: Items were adapted from the Loyalty Scale of the SUPR-Q (Sauro, 2015). CL1 used a 7-point Likert scale from 1 (strongly disagree) to 7 (completely agree); CL2 used a 7-point likelihood scale from 1 (not at all likely) to 7 (extremely likely).",
    notation = "none"
  )

Table 4. Customer Loyalty

Item	Text
CL1	I will likely visit the website in the future.
CL2	How likely are you to recommend the website to a friend or colleague?
Note: Items were adapted from the Loyalty Scale of the SUPR-Q (Sauro, 2015). CL1 used a 7-point Likert scale from 1 (strongly disagree) to 7 (completely agree); CL2 used a 7-point likelihood scale from 1 (not at all likely) to 7 (extremely likely).

We simulated customer loyalty responses as strongly influenced by both perceived ease of use and social responsibility perceptions.

# Simulate preliminary distribution using function from section 3.4
cl.scores <- f.sim.likert(n = data.n, 
                          n.items = 2,
                          mean = 4.2, 
                          sd = 1.4,
                          rounding = FALSE,
                          scale.max = 7,
                          scale.min = 1, 
                          corr = 0.6)

# Apply relationships with ease of use and social responsibility perceptions

# Define effect sizes (larger than for purchase intention)
csr.effect.size <- 0.7
eou.effect.size <- 0.6

# Apply effects to each item
for(i in 1:ncol(cl.scores)) {
  cl.scores[, i] <- cl.scores[, i] + csr.mean * csr.effect.size + eou.mean * eou.effect.size
}

# Ensure scores remain within scale limits
cl.scores <- round(pmin(pmax(cl.scores, 1), 7))

# Assign scores to dataframe
pay.df$cl1 <- cl.scores[, 1]
pay.df$cl2 <- cl.scores[, 2]

# Check distributions and relationships
par(mfrow = c(1, 2))
barplot(table(pay.df$cl1), main = "CL1 Distribution")
barplot(table(pay.df$cl2), main = "CL2 Distribution")

# Check inter-item correlations
corr.test(pay.df[, c("cl1", "cl2")])

## Call:corr.test(x = pay.df[, c("cl1", "cl2")])
## Correlation matrix 
##      cl1  cl2
## cl1 1.00 0.68
## cl2 0.68 1.00
## Sample Size 
## [1] 800
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##     cl1 cl2
## cl1   0   0
## cl2   0   0
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

# Examine scores across experimental conditions
aggregate(pay.df$cl1, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

aggregate(pay.df$cl2, by = list(condition = pay.df$condition, customer = pay.df$customer), mean)

3.4.5. Payment Method Adoption

Adoption of ControvPay was measured at the end of the survey by asking participants which payment method they would choose among the options presented. We simulated selection of ControvPay in the experimental condition, with lower probabilities for participants who provided low social responsibility scores.

# Set reference probability for adoption
prob.adoption <- rep(0.15, data.n)

# Reduce probability for participants with low social responsibility perceptions
prob.adoption[pay.df$condition == "experimental" & csr.mean <= -1] <- 0.05

# Create adoption variable
pay.df$adoption <- rbinom(n = data.n, size = 1, prob = prob.adoption)

# Set adoption to NA for control condition (ControvPay not available)
pay.df$adoption[pay.df$condition == "control"] <- NA_real_

# Check distribution
with(pay.df, table(adoption, condition))

##         condition
## adoption control experimental
##        0       0          356
##        1       0           42

4. Reliability Analyses

We assessed the internal reliability of the multi-item scales (social responsibility, ease of use, purchase intention, and loyalty) using Cronbach’s alpha.

# Social responsibility
rel.csr <- alpha(x = pay.df[, c("csr1", "csr2", "csr3", "csr4")])
summary(rel.csr)

## 
## Reliability analysis   
##  raw_alpha std.alpha G6(smc) average_r S/N  ase mean  sd median_r
##       0.82      0.82    0.78      0.54 4.7 0.01  5.4 1.1     0.54

# Ease of use
rel.eou <- alpha(x = pay.df[, c("eou1", "eou2", "eou3", "eou4")])
summary(rel.eou)

## 
## Reliability analysis   
##  raw_alpha std.alpha G6(smc) average_r S/N    ase mean  sd median_r
##       0.87      0.87    0.83      0.62 6.5 0.0077  4.2 1.4     0.61

# Purchase intention
rel.pi <- alpha(x = pay.df[, c("pi1", "pi2")])
summary(rel.pi)

## 
## Reliability analysis   
##  raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
##       0.75      0.75     0.6       0.6 3.1 0.017    5 1.1      0.6

# Customer loyalty
rel.cl <- alpha(x = pay.df[, c("cl1", "cl2")])
summary(rel.cl)

## 
## Reliability analysis   
##  raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
##       0.81      0.81    0.68      0.68 4.2 0.014  4.2 1.4     0.68

All four scales demonstrated acceptable internal reliability (Cronbach’s α > 0.70), supporting the calculation of composite scale scores as the mean of each scale’s items.

# Create composite scale scores
pay.df$social.responsibility <- rowMeans(pay.df[, c("csr1", "csr2", "csr3", "csr4")])
pay.df$ease <- rowMeans(pay.df[, c("eou1", "eou2", "eou3", "eou4")])
pay.df$purchase.intention <- rowMeans(pay.df[, c("pi1", "pi2")])
pay.df$loyalty <- rowMeans(pay.df[, c("cl1", "cl2")])

# Move individual items to separate dataframe to maintain clean primary dataset
item.names <- c("csr1", "csr2", "csr3", "csr4", 
                "eou1", "eou2", "eou3", "eou4", 
                "pi1", "pi2", 
                "cl1", "cl2")
items.df <- pay.df[, names(pay.df) %in% item.names]
pay.df <- pay.df[, !names(pay.df) %in% item.names]

# Check dataset structure
summary(pay.df)

##        id             condition        age          gender   
##  1      :  1   control     :402   Min.   :16.00   man  :399  
##  2      :  1   experimental:398   1st Qu.:32.00   woman:401  
##  3      :  1                      Median :40.50              
##  4      :  1                      Mean   :42.53              
##  5      :  1                      3rd Qu.:50.00              
##  6      :  1                      Max.   :85.00              
##  (Other):794                                                 
##           freq.customer customer     adoption      social.responsibility
##  none            :396   no :396   Min.   :0.0000   Min.   :1.750        
##  1 time          :162   yes:404   1st Qu.:0.0000   1st Qu.:4.500        
##  2-4 times       :110             Median :0.0000   Median :5.500        
##  5-9 times       : 89             Mean   :0.1055   Mean   :5.364        
##  10 times or more: 43             3rd Qu.:0.0000   3rd Qu.:6.250        
##                                   Max.   :1.0000   Max.   :7.000        
##                                   NA's   :402                           
##       ease       purchase.intention    loyalty     
##  Min.   :1.000   Min.   :1.500      Min.   :1.000  
##  1st Qu.:3.250   1st Qu.:4.000      1st Qu.:3.500  
##  Median :4.250   Median :5.000      Median :4.250  
##  Mean   :4.239   Mean   :4.967      Mean   :4.244  
##  3rd Qu.:5.250   3rd Qu.:5.500      3rd Qu.:5.500  
##  Max.   :7.000   Max.   :7.000      Max.   :7.000  
## 

5. Sample Demographics

To obtain a comprehensive sample summary, we used:

• The summary() function to examine frequency distributions of categorical variables.
• The describe() function from the psych package to obtain descriptive statistics for numeric variables.

# Descriptive statistics for numeric variables
describe(pay.df)

# Frequency distributions for categorical variables  
summary(pay.df)

##        id             condition        age          gender   
##  1      :  1   control     :402   Min.   :16.00   man  :399  
##  2      :  1   experimental:398   1st Qu.:32.00   woman:401  
##  3      :  1                      Median :40.50              
##  4      :  1                      Mean   :42.53              
##  5      :  1                      3rd Qu.:50.00              
##  6      :  1                      Max.   :85.00              
##  (Other):794                                                 
##           freq.customer customer     adoption      social.responsibility
##  none            :396   no :396   Min.   :0.0000   Min.   :1.750        
##  1 time          :162   yes:404   1st Qu.:0.0000   1st Qu.:4.500        
##  2-4 times       :110             Median :0.0000   Median :5.500        
##  5-9 times       : 89             Mean   :0.1055   Mean   :5.364        
##  10 times or more: 43             3rd Qu.:0.0000   3rd Qu.:6.250        
##                                   Max.   :1.0000   Max.   :7.000        
##                                   NA's   :402                           
##       ease       purchase.intention    loyalty     
##  Min.   :1.000   Min.   :1.500      Min.   :1.000  
##  1st Qu.:3.250   1st Qu.:4.000      1st Qu.:3.500  
##  Median :4.250   Median :5.000      Median :4.250  
##  Mean   :4.239   Mean   :4.967      Mean   :4.244  
##  3rd Qu.:5.250   3rd Qu.:5.500      3rd Qu.:5.500  
##  Max.   :7.000   Max.   :7.000      Max.   :7.000  
## 

The sample characteristics indicated that:

• Gender categories were evenly distributed, with approximately 50% men and 50% women.
• The average participant age was 43 years (SD = 14), ranging from 16 to 85 years.
• Around half of the sample consisted of the retailer’s customers (having made at least one purchase in the previous year).
• Participants were evenly assigned to control and experimental conditions, with each group comprising 50% of the sample.

However, to ensure the validity of comparisons across customer and experimental groups, we needed to check for demographic differences between these groups.

6. Checks for Balanced Comparisons

All comparisons in this study rely on groups defined by experimental assignment and customer status. It was therefore essential to examine whether demographic differences could explain potential effects between these groups.

6.1. Balance in Customer Groups

We assessed whether the retailer’s customers and non-customers differed on age and gender using descriptive statistics, visualizations, and statistical tests (independent t-tests for age, chi-square tests for gender).

# Descriptive statistics by customer status
describeBy(pay.df[, c("age", "gender")], group = pay.df$customer)

## 
##  Descriptive statistics by group 
## group: no
##        vars   n  mean    sd median trimmed   mad min max range skew kurtosis
## age       1 396 42.49 14.22     40   41.08 13.34  16  85    69 0.94     0.92
## gender    2 396  1.49  0.50      1    1.49  0.00   1   2     1 0.04    -2.00
##          se
## age    0.71
## gender 0.03
## ------------------------------------------------------------ 
## group: yes
##        vars   n  mean    sd median trimmed   mad min max range  skew kurtosis
## age       1 404 42.58 14.61     41   41.35 13.34  16  85    69  0.76     0.36
## gender    2 404  1.51  0.50      2    1.52  0.00   1   2     1 -0.05    -2.00
##          se
## age    0.73
## gender 0.02

# Statistical test for age differences
t.test(age ~ customer, data = pay.df)

## 
##  Welch Two Sample t-test
## 
## data:  age by customer
## t = -0.087598, df = 797.96, p-value = 0.9302
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -2.090580  1.911962
## sample estimates:
##  mean in group no mean in group yes 
##          42.48990          42.57921

cohens_d(age ~ customer, data = pay.df)

# Statistical test for gender distribution
table.gender.customer <- with(pay.df, table(gender, customer))
chisq.test(table.gender.customer)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table.gender.customer
## X-squared = 0.31923, df = 1, p-value = 0.5721

# Visualizations
par(mfrow = c(1, 2))
boxplot(age ~ customer, data = pay.df, main = "Age by Customer Status")
barplot(prop.table(table.gender.customer, margin = 2), beside = TRUE, 
        main = "Gender Distribution by Customer Status")

par(mfrow = c(1, 1))

Retailer customers and non-customers showed similar age distributions and gender composition, suggesting that any subsequent differences between these groups are unlikely to be attributable to these demographic factors.

6.2. Balance in Experimental Assignment

We assessed whether control and experimental groups differed on age, gender, or customer status using descriptive statistics, visualizations, and statistical tests (independent t-tests for age, chi-square tests for categorical variables).

# Descriptive statistics by experimental condition
describeBy(pay.df[, c("age", "gender", "customer")], group = pay.df$condition)

## 
##  Descriptive statistics by group 
## group: control
##          vars   n  mean    sd median trimmed   mad min max range skew kurtosis
## age         1 402 42.29 14.21     40   40.93 13.34  16  85    69 0.88     0.68
## gender      2 402  1.50  0.50      1    1.50  0.00   1   2     1 0.01    -2.00
## customer    3 402  1.48  0.50      1    1.47  0.00   1   2     1 0.09    -2.00
##            se
## age      0.71
## gender   0.02
## customer 0.02
## ------------------------------------------------------------ 
## group: experimental
##          vars   n  mean    sd median trimmed   mad min max range  skew kurtosis
## age         1 398 42.78 14.62     41   41.51 13.34  16  85    69  0.82     0.57
## gender      2 398  1.51  0.50      2    1.51  0.00   1   2     1 -0.02    -2.00
## customer    3 398  1.53  0.50      2    1.54  0.00   1   2     1 -0.13    -1.99
##            se
## age      0.73
## gender   0.03
## customer 0.03

# Statistical test for age differences
t.test(age ~ condition, data = pay.df)

## 
##  Welch Two Sample t-test
## 
## data:  age by condition
## t = -0.48576, df = 796.82, p-value = 0.6273
## alternative hypothesis: true difference in means between group control and group experimental is not equal to 0
## 95 percent confidence interval:
##  -2.497106  1.506381
## sample estimates:
##      mean in group control mean in group experimental 
##                   42.28856                   42.78392

cohens_d(age ~ condition, data = pay.df)

# Statistical test for gender distribution
table.gender.condition <- with(pay.df, table(gender, condition))
chisq.test(table.gender.condition)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table.gender.condition
## X-squared = 0.020101, df = 1, p-value = 0.8873

# Statistical test for customer distribution  
table.customer.condition <- with(pay.df, table(customer, condition))
chisq.test(table.customer.condition)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table.customer.condition
## X-squared = 2.2095, df = 1, p-value = 0.1372

# Visualizations
par(mfrow = c(1, 3))
boxplot(age ~ condition, data = pay.df, main = "Age by Condition")
barplot(prop.table(table.gender.condition, margin = 2), beside = TRUE, 
        main = "Gender Distribution by Condition")
barplot(prop.table(table.customer.condition, margin = 2), beside = TRUE, 
        main = "Customer Distribution by Condition")

par(mfrow = c(1, 1))

The two experimental groups were similar in terms of age, gender distribution, and customer ratio.

Since customer status is considered as a potential interaction factor, we also examined assignment balance within customer and non-customer subgroups.

6.2.1. Balance in Experimental Assignment within Customer Groups

We examined whether demographic characteristics differed across the four experimental cells (control/experimental × customer/non-customer) using descriptive statistics, visualizations, and statistical tests (factorial ANOVA for age, chi-square test for gender).

# Descriptive statistics by condition and customer groups
describeBy(pay.df[, c("age", "gender")], group = list(pay.df$condition, pay.df$customer))

## 
##  Descriptive statistics by group 
## : control
## : no
##        vars   n  mean    sd median trimmed   mad min max range skew kurtosis
## age       1 210 41.64 13.69     39   40.40 13.34  16  85    69 0.89     0.91
## gender    2 210  1.48  0.50      1    1.48  0.00   1   2     1 0.08    -2.00
##          se
## age    0.94
## gender 0.03
## ------------------------------------------------------------ 
## : experimental
## : no
##        vars   n  mean    sd median trimmed   mad min max range skew kurtosis
## age       1 186 43.45 14.78   41.0   41.91 11.86  16  85    69 0.97     0.79
## gender    2 186  1.50  0.50    1.5    1.50  0.74   1   2     1 0.00    -2.01
##          se
## age    1.08
## gender 0.04
## ------------------------------------------------------------ 
## : control
## : yes
##        vars   n  mean    sd median trimmed   mad min max range  skew kurtosis
## age       1 192 42.99 14.77   41.5   41.58 14.08  19  85    66  0.84     0.39
## gender    2 192  1.52  0.50    2.0    1.52  0.00   1   2     1 -0.06    -2.01
##          se
## age    1.07
## gender 0.04
## ------------------------------------------------------------ 
## : experimental
## : yes
##        vars   n  mean   sd median trimmed   mad min max range  skew kurtosis
## age       1 212 42.20 14.5     41   41.09 13.34  16  85    69  0.67     0.28
## gender    2 212  1.51  0.5      2    1.51  0.00   1   2     1 -0.04    -2.01
##          se
## age    1.00
## gender 0.03

# Factorial ANOVA for age differences
age_anova <- aov(age ~ condition * customer, data = pay.df)
summary(age_anova)

##                     Df Sum Sq Mean Sq F value Pr(>F)
## condition            1     49    49.1   0.236  0.627
## customer             1      1     0.8   0.004  0.951
## condition:customer   1    336   335.7   1.614  0.204
## Residuals          796 165583   208.0

age.eta.squared <- eta_squared(age_anova, partial = TRUE) 
age.eta.squared

# Create combined grouping variable for chi-square test
pay.df$condition.customer <- rep(NA_character_, nrow(pay.df))
pay.df$condition.customer[pay.df$condition == "experimental" & pay.df$customer == "yes"] <- "exp.customers"
pay.df$condition.customer[pay.df$condition == "experimental" & pay.df$customer == "no"] <- "exp.non.customers"
pay.df$condition.customer[pay.df$condition == "control" & pay.df$customer == "yes"] <- "ct.customers"
pay.df$condition.customer[pay.df$condition == "control" & pay.df$customer == "no"] <- "ct.non.customers"
pay.df$condition.customer <- factor(pay.df$condition.customer)

# Chi-square test for gender distribution
table.gender.condition.customer <- with(pay.df, table(gender, condition.customer))
chisq.test(table.gender.condition.customer)

## 
##  Pearson's Chi-squared test
## 
## data:  table.gender.condition.customer
## X-squared = 0.56274, df = 3, p-value = 0.9049

# Visualizations for the four groups
par(mar = c(8, 4.1, 4.1, 2.1))  # Adjust margins for x-axis labels

# Age distribution
boxplot(age ~ condition.customer, data = pay.df,
        main = "Age Distribution by Condition × Customer",
        las = 2)  # Rotate x-axis labels

# Gender proportions
barplot(prop.table(table.gender.condition.customer, margin = 2), beside = TRUE,
        main = "Gender Distribution by Condition × Customer",
        legend.text = TRUE)

par(mar = c(5.1, 4.1, 4.1, 2.1))  # Reset margins

The factorial ANOVA revealed no significant interaction between condition and customer status for age, and the chi-square test showed no significant association between the combined condition-customer groups and gender distribution. This indicates successful balance across all four experimental cells, supporting the internal validity of subsequent comparisons.

Given the minimal demographic differences across our experimental factors, we proceeded with subsequent analyses without including these variables as covariates.

7. Preliminary Exploration of Data Correlations

To gain a more holistic understanding of variable interrelationships, we examined correlations among demographic, experimental, and outcome variables using correlation matrices, scatterplot matrices, and correlation plots.

# Create dataframe with variables for correlation analysis
cor.df <- pay.df[, c("age", "gender", "condition", "customer", 
                     "social.responsibility", "ease", "purchase.intention", 
                     "loyalty", "adoption")]  # Analyze adoption separately due to missing data pattern

# Convert categorical variables to numeric for correlation analysis
cor.df$gender <- as.numeric(ifelse(cor.df$gender == "woman", 0, 1))
cor.df$condition <- as.numeric(ifelse(cor.df$condition == "control", 0, 1))
cor.df$customer <- as.numeric(ifelse(cor.df$customer == "no", 0, 1))

# Correlation matrix using pairwise deletion
corr.test(cor.df, use = "pairwise")

## Warning in cor(x, use = use, method = method): La desviación estándar es cero

## Call:corr.test(x = cor.df, use = "pairwise")
## Correlation matrix 
##                         age gender condition customer social.responsibility
## age                    1.00   0.02      0.02     0.00                 -0.04
## gender                 0.02   1.00     -0.01    -0.02                 -0.05
## condition              0.02  -0.01      1.00     0.06                 -0.06
## customer               0.00  -0.02      0.06     1.00                 -0.04
## social.responsibility -0.04  -0.05     -0.06    -0.04                  1.00
## ease                  -0.37  -0.03      0.00    -0.01                  0.04
## purchase.intention    -0.08  -0.04     -0.01     0.02                  0.11
## loyalty               -0.18   0.00     -0.04    -0.04                  0.44
## adoption               0.05  -0.03        NA     0.04                  0.06
##                        ease purchase.intention loyalty adoption
## age                   -0.37              -0.08   -0.18     0.05
## gender                -0.03              -0.04    0.00    -0.03
## condition              0.00              -0.01   -0.04       NA
## customer              -0.01               0.02   -0.04     0.04
## social.responsibility  0.04               0.11    0.44     0.06
## ease                   1.00               0.23    0.38    -0.04
## purchase.intention     0.23               1.00    0.11     0.04
## loyalty                0.38               0.11    1.00    -0.02
## adoption              -0.04               0.04   -0.02     1.00
## Sample Size 
##                       age gender condition customer social.responsibility ease
## age                   800    800       800      800                   800  800
## gender                800    800       800      800                   800  800
## condition             800    800       800      800                   800  800
## customer              800    800       800      800                   800  800
## social.responsibility 800    800       800      800                   800  800
## ease                  800    800       800      800                   800  800
## purchase.intention    800    800       800      800                   800  800
## loyalty               800    800       800      800                   800  800
## adoption              398    398       398      398                   398  398
##                       purchase.intention loyalty adoption
## age                                  800     800      398
## gender                               800     800      398
## condition                            800     800      398
## customer                             800     800      398
## social.responsibility                800     800      398
## ease                                 800     800      398
## purchase.intention                   800     800      398
## loyalty                              800     800      398
## adoption                             398     398      398
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##                        age gender condition customer social.responsibility ease
## age                   0.00   1.00      1.00     1.00                  1.00 0.00
## gender                0.50   0.00      1.00     1.00                  1.00 1.00
## condition             0.63   0.83      0.00     1.00                  1.00 1.00
## customer              0.93   0.53      0.12     0.00                  1.00 1.00
## social.responsibility 0.26   0.16      0.07     0.29                  0.00 1.00
## ease                  0.00   0.47      0.94     0.69                  0.23 0.00
## purchase.intention    0.03   0.29      0.81     0.66                  0.00 0.00
## loyalty               0.00   0.91      0.20     0.27                  0.00 0.00
## adoption              0.29   0.56        NA     0.39                  0.23 0.46
##                       purchase.intention loyalty adoption
## age                                 0.77    0.00        1
## gender                              1.00    1.00        1
## condition                           1.00    1.00       NA
## customer                            1.00    1.00        1
## social.responsibility               0.07    0.00        1
## ease                                0.00    0.00        1
## purchase.intention                  0.00    0.07        1
## loyalty                             0.00    0.00        1
## adoption                            0.48    0.65        0
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

# Scatterplot matrix with continuous variables only
scatterplotMatrix(~ age + social.responsibility + ease + purchase.intention + loyalty, 
                  data = pay.df, use = "complete.obs")

# Correlation plot using pairwise correlations
corrplot.mixed(corr = cor(cor.df, use = "pairwise.complete.obs"), upper = "ellipse")

## Warning in cor(cor.df, use = "pairwise.complete.obs"): La desviación estándar
## es cero

The correlation analysis revealed several key relationships:

• Age showed moderate negative associations with ease of use, purchase intention, and loyalty. However, since experimental conditions and customer groups demonstrated balanced age distributions, we proceeded with primary analyses without including age as a covariate.
• Experimental exposure to ControvPay showed minimal correlations with dependent variables in bivariate analyses. We therefore proceeded to more complex models incorporating interaction effects with customer status to detect potential conditional effects (Section 9).
• Social responsibility and ease perceptions were positively associated with higher purchase intention and loyalty. These relationships align with theoretical expectations and support the specification of mediation models to examine underlying mechanisms (Section 11).

8. Expected Adoption of ControvPay

To evaluate whether ControvPay adoption differs between the retailer’s customers and non-customers, we conducted a chi-square test and examined frequency distributions.

# Frequency table
adoption.customer.table <- with(pay.df, table(adoption, customer))
adoption.customer.table

##         customer
## adoption  no yes
##        0 169 187
##        1  17  25

# Proportion table
prop.table(adoption.customer.table, margin = 2)

##         customer
## adoption         no        yes
##        0 0.90860215 0.88207547
##        1 0.09139785 0.11792453

# Chi-square test
chisq.test(adoption.customer.table)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  adoption.customer.table
## X-squared = 0.48429, df = 1, p-value = 0.4865

The results indicate no significant difference in adoption rates between current customers and non-customers. We therefore estimated overall ControvPay adoption across the entire sample.

# Overall frequency table
adoption.table <- with(pay.df, table(adoption))

# Proportion table
prop.table(adoption.table)

## adoption
##         0         1 
## 0.8944724 0.1055276

# Confidence intervals for adoption rate
adoption.prop <- binom.test(x = adoption.table["1"],
                            n = sum(adoption.table))

# Prepare results dataframe
adoption.df <- data.frame(
  category = c("did not adopt ControvPay", "adopted ControvPay"),
  percentage = as.numeric(prop.table(adoption.table)) * 100,
  lower.ci = NA_real_,
  upper.ci = NA_real_,
  count = as.numeric(adoption.table)
)
adoption.df$lower.ci[adoption.df$category == "adopted ControvPay"] <- adoption.prop$conf.int[1] * 100
adoption.df$upper.ci[adoption.df$category == "adopted ControvPay"] <- adoption.prop$conf.int[2] * 100

# Display results
adoption.df

# Plotting the results in a bar graph
ggplot(adoption.df[adoption.df$category == "adopted ControvPay", ], aes(x = category, y = percentage)) +
  
  # Type of graph
  geom_col(width = 0.7, fill = "grey80") +
  
  # Error bars
  geom_errorbar(
    aes(ymin = lower.ci, ymax = upper.ci),
        width = 0.2) +
  
  # Graph title
  labs(
    title = "Adoption of ControvPay",
    y = "% of Customers",
    x = " " ) +
  
  # Data labels
  geom_text(
    aes(y = upper.ci, label = paste0(round(percentage, 0), "%")),         # label with percentage rounded
    position = position_dodge(width = 0.8),   # match the dodge of bars
    vjust = -1,                             # slightly above the bar
    size = 5
  ) +
  
  # Y axis limits para que quepan las etiquetas
  ylim(0, max(adoption.df$upper.ci, na.rm = TRUE) * 1.2) +
  
  # Formatting
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"), 
           axis.text.x = element_blank() 
  )

The results indicate that 11% of users would choose ControvPay as their preferred payment option if available.

While this adoption rate supports including ControvPay to accommodate diverse customer preferences, subsequent analyses will examine its potential impact on overall payment process perceptions.

9. Experimental Effects Analysis

To evaluate the causal impact of exposing participants to ControvPay versus traditional payment methods only, we employed 2×2 factorial Analysis of Variance (ANOVA). This design enabled testing both the main effect of experimental condition and its potential interaction with customer status (current customers vs. non-customers).

The analysis addresses two key research questions:

1. Does the presence of ControvPay influence consumer perceptions and intentions?
2. Do these effects differ between existing customers and non-customers?

We began by examining the normality assumption for ANOVA tests.

9.1. Checking the Normality Assumption

To evaluate the normality assumption for ANOVA, we examined histograms, normal Q-Q plots, and calculated skewness and kurtosis indices.

# Skewness and kurtosis indices
normality.indexes.df <- describe(pay.df[, c("social.responsibility", "ease", "purchase.intention", "loyalty")])
normality.indexes.df <- normality.indexes.df[, c("skew", "kurtosis")]
normality.indexes.df 

# Prepare layout for 4 plots
par(mfrow = c(2, 2))

# Histograms
hist(pay.df$social.responsibility, main = "Social Responsibility")
hist(pay.df$ease, main = "Ease of Use")
hist(pay.df$purchase.intention, main = "Purchase Intention")
hist(pay.df$loyalty, main = "Loyalty")

# Q-Q plots
qqnorm(pay.df$social.responsibility, main = "Social Responsibility Q-Q Plot")
qqline(pay.df$social.responsibility)

qqnorm(pay.df$ease, main = "Ease of Use Q-Q Plot")
qqline(pay.df$ease)

qqnorm(pay.df$purchase.intention, main = "Purchase Intention Q-Q Plot")
qqline(pay.df$purchase.intention)

qqnorm(pay.df$loyalty, main = "Loyalty Q-Q Plot")
qqline(pay.df$loyalty)

The distributions of ease of use perceptions, purchase intention and loyalty closely approximate theoretical normal distributions. The histograms for social responsibility perceptions suggest a ceiling effect, resulting in some skewness to the left. However, this level of skewness appears tolerable, as indicated by skewness indices within the ±1 range and Q-Q plots showing reasonable approximation to normality. We proceeded with ANOVA models for all variables, as ANOVA is robust to moderate normality violations with large sample sizes.

9.2. Factorial ANOVA Models

To compare conditions, we employed a 2 (condition: experimental vs. control) × 2 (customer status: customer vs. non-customer) factorial ANOVA for each dependent variable. For this analysis, we utilized a reusable function from my analytical toolkit, specified below.

# Function for obtaining factorial ANOVA results for multiple DVs across 2 factors

# Function name: f.factorial.anova.2f
#
# Input parameters:
#   my.data: the dataframe
#   my.dv.list: vector of dependent variable names
#   my.factor1: first factor variable
#   my.factor2: second factor variable
#
# Output: list containing:
#   anova: ANOVA results for main effects and interactions, with effect sizes (eta squared)
#   descriptives: means, standard deviations, and 95% CI for all effects
#   tukey: post-hoc comparisons with Cohen's d effect sizes

f.factorial.anova.2f <- function (my.data, my.dv.list, my.factor1, my.factor2) {
  
  # Initialize lists for results storage
  anova.results <- list()    # ANOVA results
  descriptive.results <- list()  # Descriptive statistics
  cohen.d.results <- list()  # Cohen's d results
  tukey.results <- list()    # Tukey post-hoc results
  a <- d <- c <- t <- 1  # List indices for the lists above

  # Obtain factor levels
  f1.levels <- unique (my.data[[my.factor1]])
  f2.levels <- unique (my.data[[my.factor2]])
  
  # Specify function for Cohen d calculation
    f.cohen.d.calculation <- function (groups.data, groups.matrix, dv.name, effect.name, level.names) {
    
      # Create a list to store the results
      fd.results.list <- list()
      fd <- 1 # index for the above list
    
      # Open a loop for each comparison
      for (i.comparison in 1:nrow(groups.matrix)){
      
        # Obtain the group indexes to compare
        index.g1 <- groups.matrix [i.comparison, 1]
        index.g2 <- groups.matrix [i.comparison, 2]
      
        # Obtain the names of groups compared
        name.g1 <- groups.data [index.g1, level.names]
        name.g2 <- groups.data [index.g2, level.names]
      
        # Obtain the means of groups compared
        mean.g1 <- groups.data [index.g1, "mean"]
        mean.g2 <- groups.data [index.g2, "mean"]
      
        # Obtain the sd of groups compared
        sd.g1 <- groups.data [index.g1, "sd"]
        sd.g2 <- groups.data [index.g2, "sd"]
      
        # Obtain the n of groups compared
        n.g1 <- groups.data [index.g1, "n"]
        n.g2 <- groups.data [index.g2, "n"]     
      
        # Calculate the pooled standard deviation
        pooled.sd <- sqrt(((n.g1 - 1) * sd.g1^2 + (n.g2 - 1) * sd.g2^2) / (n.g1 + n.g2 - 2))
      
        # Calculate Cohen's d
        cohen.d <- (mean.g1 - mean.g2)/ pooled.sd
      
        # Save the results for each comparison
        comparison.cohen.d.result <- data.frame (
          dep.variable = dv.name,
          effect = effect.name,
          group1 = name.g1,
          group2 = name.g2,
          d = cohen.d, 
          sd.pooled = pooled.sd
        )
        fd.results.list [[fd]] <- comparison.cohen.d.result
        fd <- fd + 1
      
      } # Close the loop within the f.cohen.d.calculation   
      
      d.function.results <- do.call(rbind, fd.results.list)
      
    # Return the result for the effect
    return (d.function.results)
   
    } # Close the f.cohen.d.calculation function
      
  # Initiate a loop for each dv
  for (dv in my.dv.list) {
    
    
#--- 1. Metrics for the ANOVA models
    
    # Build a dynamic formula
    anova.formula <- as.formula (paste0(dv, "~", my.factor1, "*", my.factor2))

    # Fit the ANOVA model
    anova.model <- aov(anova.formula, data = my.data)
    
    # Extract the model in a dataframe
    anova.table <- anova(anova.model)
    
    # Bind the anova table with a dataframe that contains the dv names
    anova.result.df <- data.frame (dep.variable = rep(dv, nrow(anova.table)), 
                                    effect = rownames(anova.table))
    anova.result.df <- bind_cols (anova.result.df, anova.table)
    
    # Add the eta squared (Proportional Reduction in Error)
    ss.total <- sum(anova.table[ , "Sum Sq"]) # sum of error
    anova.result.df$eta.squared <- anova.result.df[["Sum Sq"]] / ss.total
    
    # Save the results
    anova.results [[a]] <- anova.result.df
    a <- a + 1
    
 
    
#--- 2. Posthoc Comparison with Tukey
    
    # Obtain the Tukey results
    tukey <- TukeyHSD(anova.model)
    
    # Bind the comparisons for main effects and interaction effects in a single dataframe
    for (name.effect in names(tukey)) {
      
      # Extract the comparisons for the effect
      tukey.df <- as.data.frame(tukey[[name.effect]])
      tukey.df$dep.variable <- dv
      tukey.df$effect <- name.effect
      tukey.df$comparison <- rownames(tukey.df)
      tukey.results [[t]] <- tukey.df
      t <- t + 1
    } # Close the loop for the Tukey effect       
    
    
    
#--- 3. Descriptive Statistics
    
    # -- 3.1. Descriptives for factor 1
    
    # Open a loop for factor 1
    for (f1.level in f1.levels){
      
      # Create a vector with the scores of each level in factor 1
      current.vector <- my.data[my.data[[my.factor1]] == f1.level , dv]
      
      # Clean missind data
      current.vector.c <- current.vector[!is.na(current.vector)]
      
      # Obtain the mean, standard deviations, standard errors, n, and 95% CI
      mean.level <- mean(current.vector.c)
      sd.level <- sd(current.vector.c)
      n.level <- length(current.vector.c)
      se.level <- sd.level / sqrt(n.level)
      upper.ci.level <- mean.level + 1.96 * se.level
      lower.ci.level <- mean.level - 1.96 * se.level
      
      # Save the results
      level.descriptive.results <- data.frame (
        dep.variable = dv, 
        factor1.levels = f1.level,
        factor2.levels = paste0 ("all"),
        mean = mean.level,
        sd = sd.level,
        upper.ci = upper.ci.level,
        lower.ci = lower.ci.level,
        se = se.level,
        n = n.level, 
        level.name = f1.level)
      
      
      descriptive.results[[d]] <- level.descriptive.results
      d <- d + 1
    
    } # Close the loop for factor 1
     
    
    # -- 3.2. Descriptives for factor 2
      
    # Open a loop for factor 2
    for (f2.level in f2.levels){
      
      # Create a vector with the scores of each level in factor 2
      current.vector <- my.data[my.data[[my.factor2]] == f2.level , dv]
      
      # Clean missind data
      current.vector.c <- current.vector[!is.na(current.vector)]
      
      # Obtain the mean, standard deviations, standard errors, n, and 95% CI
      mean.level <- mean(current.vector.c)
      sd.level <- sd(current.vector.c)
      n.level <- length(current.vector.c)
      se.level <- sd.level / sqrt(n.level)
      upper.ci.level <- mean.level + 1.96 * se.level
      lower.ci.level <- mean.level - 1.96 * se.level
      
      # Save the results
      level.descriptive.results <- data.frame (
        dep.variable = dv, 
        factor1.levels = paste0 ("all"),
        factor2.levels = f2.level,
        mean = mean.level,
        sd = sd.level,
        upper.ci = upper.ci.level,
        lower.ci = lower.ci.level,
        se = se.level,
        n = n.level,
        level.name = f2.level
      )
      
      descriptive.results[[d]] <- level.descriptive.results
      d <- d + 1
      
      } # Close the loop for factor 2
    
    
    # -- 3.3. Descriptives for the interaction
    
      # Open loops for the interaction

      for (f1.level in f1.levels){
        for (f2.level in f2.levels){
      
        # Create a vector with the scores of the interaction level
        current.vector <- my.data [my.data[[my.factor1]] == f1.level &  my.data[[my.factor2]] == f2.level, dv ]
      
        # Clean missind data
        current.vector.c <- current.vector[!is.na(current.vector)]
      
        # Obtain the mean, standard deviations, standard errors, n, and 95% CI
        mean.level <- mean(current.vector.c)
        sd.level <- sd(current.vector.c)
        n.level <- length(current.vector.c)
        se.level <- sd.level / sqrt(n.level)
        upper.ci.level <- mean.level + 1.96 * se.level
        lower.ci.level <- mean.level - 1.96 * se.level
      
        # Save the results
        level.descriptive.results <- data.frame (
        dep.variable = dv, 
        factor1.levels = f1.level,
        factor2.levels = f2.level,
        mean = mean.level,
        sd = sd.level,
        upper.ci = upper.ci.level,
        lower.ci = lower.ci.level,
        se = se.level,
        n = n.level,
        level.name = paste0(f1.level, ":", f2.level)
      )
      
      descriptive.results[[d]] <- level.descriptive.results
      d <- d + 1
    
    } # Close the loop of factor 1 in the interaction loops
      }  # Close the loop of factor 2 in the interaction loop
    
    
# --- 4. Cohen's d Calculations
    
    # -- 4.1. Cohen's d for factor 1 comparisons
    
    # Extract the descriptives for factor 1
    data.groups.to.compare <- do.call(rbind, descriptive.results)
    data.groups.to.compare <- data.groups.to.compare[data.groups.to.compare$dep.variable == dv & 
                                   data.groups.to.compare$factor2.levels == "all", ]
    
    # Obtain the number of groups
    n.groups <- nrow(data.groups.to.compare)
    
    # Obtain a matrix with the group combination indexes
    comparison.matrix <- as.data.frame(t(combn(1:n.groups, 2)))
    
    # Apply the f.cohen.d.calculation function to calculate the d results for factor 1
    factor1.d.results <- f.cohen.d.calculation(
      groups.data = data.groups.to.compare,
      groups.matrix = comparison.matrix,
      dv.name = dv,
      effect.name = my.factor1, 
      level.names = "level.name"
    )
    
    # Save the d results for the factor 1  groups
    cohen.d.results[[c]] <- factor1.d.results
    c <- c + 1
    
    
    # -- 4.2. Cohen's d for factor 2 comparisons
    
    # Extract the descriptives for factor 2
    data.groups.to.compare <- do.call(rbind, descriptive.results)
    data.groups.to.compare <- data.groups.to.compare[data.groups.to.compare$dep.variable == dv & 
                                   data.groups.to.compare$factor1.levels == "all", ]
    
    # Obtain the number of groups
    n.groups <- nrow(data.groups.to.compare)
    
    # Obtain a matrix with the group combination indexes
    comparison.matrix <- as.data.frame(t(combn(1:n.groups, 2)))
    
    # Apply the f.cohen.d.calculation function to calculate the d results for factor 2
    factor2.d.results <- f.cohen.d.calculation(
      groups.data = data.groups.to.compare,
      groups.matrix = comparison.matrix,
      dv.name = dv,
      effect.name = my.factor2,
      level.names = "level.name"
    )    
    
    # Save the d results for the factor 2  groups
    cohen.d.results[[c]] <- factor2.d.results
    c <- c + 1  
    

    # -- 4.3. Cohen's d for the interaction comparisons
    
    # Extract the descriptives for the interaction
    data.groups.to.compare <- do.call(rbind, descriptive.results)
    data.groups.to.compare <- data.groups.to.compare[data.groups.to.compare$dep.variable == dv & 
                                  data.groups.to.compare$factor1.levels != "all" &
                                  data.groups.to.compare$factor2.levels != "all", 
                                  ]
    
    # Obtain the number of groups
    n.groups <- nrow(data.groups.to.compare)
    
    # Obtain a matrix with the group combination indexes
    comparison.matrix <- as.data.frame(t(combn(1:n.groups, 2)))
    
    # Apply the f.cohen.d.calculation function to calculate the d results for the interaction
    interaction.d.results <- f.cohen.d.calculation(
      groups.data = data.groups.to.compare,
      groups.matrix = comparison.matrix,
      dv.name = dv,
      effect.name = paste0(my.factor1, ":", my.factor2),
      level.names = "level.name"
    )    
    
    # Save the d results for the factor 2  groups
    cohen.d.results[[c]] <- interaction.d.results
    c <- c + 1  
    

    
} # close the dv loop
  

    
#--- 5. Save all the results

  
  #--- ANOVA Results
  
  # Group the anova results in a dataframe
  anova.results.df <- do.call(rbind, anova.results)
  rownames(anova.results.df) <- NULL
  
  # Round the data columns
  numeric.cols <- sapply(anova.results.df, is.numeric)
  anova.results.df[numeric.cols] <- round(anova.results.df[numeric.cols], 3)
  
  
  #--- Descriptive Results
  
  # Group the descriptives for all levels
  descriptives.df <- do.call(rbind, descriptive.results)
  
  # Round the columns
  numeric.cols <- sapply(descriptives.df, is.numeric)
  descriptives.df[numeric.cols] <- round(descriptives.df[numeric.cols], 2)
  
  # Rename the factor columns
  names (descriptives.df)[names(descriptives.df) == "factor1.levels"] <- my.factor1
  names (descriptives.df)[names(descriptives.df) == "factor2.levels"] <- my.factor2
  
  
  #--- Cohen's d Results
  
  # Group the d results for all comparisons
  cohen.d.df <- do.call(rbind, cohen.d.results)
  
  # Add a "comparison" variable" to use as a key for the Tukey results, duplicating all comparisons with the altered order
  cohen.d.df1 <- cohen.d.df
  cohen.d.df1$comparison <- paste0 (cohen.d.df1$group1, "-", cohen.d.df1$group2)
  cohen.d.df2 <- cohen.d.df
  cohen.d.df2$comparison <- paste0 (cohen.d.df2$group2, "-", cohen.d.df2$group1)
  cohen.d.df2$d <- cohen.d.df2$d *(-1)
  doubled.cohen.d.df <- bind_rows(cohen.d.df1, cohen.d.df2)
  
  
  #--- Tukey Results
  
  # Group all the Tukey results
  tukey.results.df <- do.call(rbind, tukey.results)
  rownames(tukey.results.df) <- NULL
  
  # Reorder the columns
  index.tukey <- tukey.results.df [ , c("dep.variable", "effect", "comparison")]
  tukey.results.df[ , c("dep.variable", "effect", "comparison")] <- NULL
  tukey.results.df <- bind_cols(index.tukey, tukey.results.df)
  
  # Add the Cohen's d results
  tukey.results.df <- merge(tukey.results.df, doubled.cohen.d.df, by = c("dep.variable", "effect", "comparison"), all.x = TRUE)
  tukey.results.df$group1 <- NULL
  tukey.results.df$group2 <- NULL

  #--- List with all the results
  
  results <- list(anova = anova.results.df, 
                  descriptives = descriptives.df,
                  cohen.d = cohen.d.df,
                  tukey = tukey.results.df)
  
  # Return the results
  return(results)

}

We applied the function to our dependent variables using experimental condition and customer status as factors.

# Vector of dependent variable names
dv.variables <- c("social.responsibility", "ease", "purchase.intention", "loyalty")

# Application of ANOVA function
pay.results <- f.factorial.anova.2f(
  my.data = pay.df,
  my.dv.list = dv.variables,
  my.factor1 = "condition",
  my.factor2 = "customer"
)

# Extract ANOVA results
pay.anova.df <- pay.results$anova
pay.anova.df

# Extract descriptive statistics
pay.descriptives.df <- pay.results$descriptives
pay.descriptives.df

# Extract post-hoc comparisons
pay.tukey.df <- pay.results$tukey
pay.tukey.df

Although no effects reached significance at p < .05, the interaction between experimental condition and customer status approached this threshold. Follow-up Tukey tests revealed that retailer customers exposed to ControvPay showed lower social responsibility perceptions compared to those not exposed (d = -0.28, Tukey p = .04).

For other dependent variables, all differences were statistically non-significant. However, descriptive patterns indicated that loyalty followed a similar trend to social responsibility, with customers showing lower loyalty when exposed to ControvPay (d = -0.12, Tukey p = .594).

We visualized these interaction effects using bar graphs:

# Select interaction descriptives for graphing
pay.descriptives.int.df <- pay.descriptives.df [!(pay.descriptives.df$condition == "all") & !(pay.descriptives.df$customer =="all"), ]

# Specify the order of the DVs
pay.descriptives.int.df$dep.variable <- factor(pay.descriptives.int.df$dep.variable, levels = c("ease", "social.responsibility", "purchase.intention", "loyalty"))

# Custom titles for faceted graphs
custom_titles <- c(
  social.responsibility   = "Social Responsibility Perception",
  ease     = "Ease Perception",
  purchase.intention = "Purchase Intention",
  loyalty   = "Loyalty"
)

# Bar graph visualization
ggplot(pay.descriptives.int.df, aes(x = condition, y = mean, fill = customer)) + 
  
  # Type of graph
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +  
  
  # Error bars
  geom_errorbar(
    aes(ymin = lower.ci, ymax = upper.ci),
    position = position_dodge(width = 0.8),
    width = 0.2
  ) + 
  
  # Facets
  facet_wrap(~ dep.variable, 
             scales = "fixed", 
             axes = "all", #
             labeller = labeller(dep.variable = custom_titles)) + 
  
  # Graph title, axis titles, and fill title
  labs(
    y = "Mean Scores",
    x = "Experimental Condition",
    fill = "Customer",
    title = "Differences Across Experimental Conditions"
  ) + 
  
  # Labels in x axis
  scale_x_discrete(labels = c ("No ControvPay", "ControvPay")) + 
  
  # Data labels
  geom_text(
    aes(y = upper.ci, label = paste0(round(mean, 2))), # use upper ci as reference position
    position = position_dodge(width = 0.8),   # match the dodge of bars
    vjust = -0.5,                             # slightly above the bar
    size = 2.5                                # text size
  ) + 
  
  # Limits for the y axis
  coord_cartesian(ylim = c(3, 6.2)) +  
  
  # Formatting
  theme_minimal(base_size = 12) +
    theme(
    axis.title.x = element_blank(),
    legend.position = "right",
    legend.justification = c("right", "bottom"),
    legend.title = element_text(size = 10, face = "bold"),
    legend.text  = element_text(size = 10),
    panel.spacing = unit(1, "lines"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    strip.text = element_text(face = "bold"),
    axis.text.x = element_text(angle = 0, vjust = 0, hjust = 0.5),  
    plot.margin = margin(10, 5, 10, 5)
  ) 

The graphs reveal small differences in social responsibility perception and loyalty among customers exposed to ControvPay. Although these effects are modest, they warrant further investigation due to their potential business implications regarding customer retention. We proceeded with Bayesian analyses to evaluate these effects more sensitively and complement our frequentist ANOVA results.

10. Bayesian Sensitivity Analysis

To conduct Bayesian ANOVA for social responsibility and loyalty perceptions, we employed a custom function from my analytical toolkit.

# Function for Bayesian ANOVA with multiple priors

# Function name: f.bayes.anova
#
# Input parameters:
#   my.data: dataframe containing the variables
#   my.dv.list: vector of dependent variable names
#   my.factor.list: vector of factor variable names  
#   my.priors: vector of prior scales to test
#
# Output: list containing:
#   individual models for each DV and prior combination
#   summary: consolidated results dataframe

f.bayes.anova <- function (my.data, my.dv.list, my.factor.list, my.priors){
  
  # Initialize results storage
  bayes.results <- list()  
  summary.bayes.results <- list()
  p <- 1 # index for summary.bayes.results
  
  # Construct formula from factor list
  factor.formula = paste0 ("~ ")
    for (i.factor in seq_along(my.factor.list)){
      if (i.factor != length(my.factor.list)) {
           factor.formula <-  paste0 (factor.formula, my.factor.list[i.factor], " * ")      
      } # Close the "if" function
           else {
             factor.formula <-  paste0 (factor.formula, my.factor.list[i.factor])
           }  # Close the "else" function
  } # Close the loop for the factors
  
  
  # Loop through DVs and priors
  for (dv in my.dv.list) {
    
    for(prior in my.priors) {
      
      # Fit Bayesian ANOVA model
      bayes.model <- anovaBF(as.formula(paste(dv, factor.formula)),
                      data = my.data,
                      rscaleFixed = prior,  
                      whichModels = 'all'
      )
      
      # Store model with descriptive name
      model.name <- paste0  (dv, ".prior.", prior)
      bayes.results[[model.name]] <- bayes.model
      
      # Create summary dataframe
      model.df <- as.data.frame (bayes.model)
      model.df$dep.variable <- dv
      model.df$prior <- prior
      model.df$effect <- row.names(model.df)
      model.df$BF <- model.df[[1]]
      model.df$error <- model.df[[2]]
      model.df <- model.df [ , c("dep.variable", "prior", "effect", "BF", "error")]
      
      # Save the summary dataframe
      summary.bayes.results [[p]] <- model.df
      p <- p  + 1 
      
    } # Close the loop for the priors
    
  } # Close the loop for the DVs
  
  # Consolidate results
  summary.df <- do.call(rbind, summary.bayes.results)
  row.names(summary.df) <- NULL
  bayes.results[["summary"]] <- summary.df
  
  return(bayes.results)
  
}

We applied the Bayesian ANOVA function to social responsibility and loyalty variables, using experimental condition and customer status as factors. To enhance sensitivity for detecting small effects, we tested multiple prior scales (r = 0.1, 0.5, 1), where r = 0.1 provides greater sensitivity to subtle effects.

# Specify priors for sensitivity analysis
priors <- c(0.1, 0.5, 1)

# Dependent variables for Bayesian analysis
dvs <- c("social.responsibility", "loyalty")

# Experimental factors
factors <- c("condition", "customer")

# Apply Bayesian ANOVA function
bayes.anova.results <- f.bayes.anova(my.data = pay.df, 
                                    my.dv.list = dvs, 
                                    my.factor.list = factors, 
                                    my.priors = priors)

# Display results summary
summary.bayes <- bayes.anova.results[["summary"]]
summary.bayes

For loyalty, the evidence consistently supported the absence of an interaction effect across all priors. While this conclusion is currently supported by anecdotal to moderate evidence, it suggests that ControvPay is unlikely to substantially impact customer loyalty.

For social responsibility, the results present ambiguous evidence regarding an interaction effect. With more conservative priors (rscale = 0.5 or 1.0), the evidence favored the null hypothesis of no differential effect between experimental conditions across customer profiles. However, these priors are less sensitive to detecting small effects. Conversely, when using a more sensitive prior (rscale = 0.1), we found anecdotal evidence for an interaction effect, though such sensitive priors also increase the risk of detecting statistical noise.

Given that even small decreases in social responsibility perceptions could be impactful for our large-scale retailer, we proceeded to quantify this potential effect using Monte Carlo simulation to estimate the practical significance of this ambiguous finding.

# Extract the model for social responsibility with sensitive prior
m.soc.resp.0.1 <- bayes.anova.results[["social.responsibility.prior.0.1"]]

# Simulate Monte Carlo chains
chain.bayes.social.responsibility <- posterior(m.soc.resp.0.1, 
                                               index = 7, 
                                               iterations = 10000)

# Convergence diagnostics for Monte Carlo chains
plot(chain.bayes.social.responsibility[, 1:10])

# Inspect Monte Carlo chains
summary(chain.bayes.social.responsibility)

## 
## Iterations = 1:10000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                                           Mean      SD  Naive SE Time-series SE
## mu                                     5.36697 0.03734 0.0003734      0.0003734
## condition-control                      0.05207 0.03501 0.0003501      0.0003847
## condition-experimental                -0.05207 0.03501 0.0003501      0.0003847
## customer-no                            0.02815 0.03322 0.0003322      0.0003471
## customer-yes                          -0.02815 0.03322 0.0003322      0.0003471
## condition:customer-control.&.no       -0.05212 0.03507 0.0003507      0.0004173
## condition:customer-control.&.yes       0.05212 0.03507 0.0003507      0.0004173
## condition:customer-experimental.&.no   0.05212 0.03507 0.0003507      0.0004173
## condition:customer-experimental.&.yes -0.05212 0.03507 0.0003507      0.0004173
## sig2                                   1.11056 0.05586 0.0005586      0.0005586
## g_condition                            0.17289 8.92045 0.0892045      0.0892045
## g_customer                             0.05757 0.85313 0.0085313      0.0085313
## g_condition:customer                   0.14459 3.02689 0.0302689      0.0302689
## 
## 2. Quantiles for each variable:
## 
##                                            2.5%       25%       50%       75%
## mu                                     5.293982  5.341575  5.366784  5.392235
## condition-control                     -0.013532  0.027819  0.051362  0.075065
## condition-experimental                -0.124319 -0.075065 -0.051362 -0.027819
## customer-no                           -0.034763  0.005577  0.027313  0.049926
## customer-yes                          -0.094969 -0.049926 -0.027313 -0.005577
## condition:customer-control.&.no       -0.125005 -0.075201 -0.050370 -0.027159
## condition:customer-control.&.yes      -0.011518  0.027159  0.050370  0.075201
## condition:customer-experimental.&.no  -0.011518  0.027159  0.050370  0.075201
## condition:customer-experimental.&.yes -0.125005 -0.075201 -0.050370 -0.027159
## sig2                                   1.007203  1.072083  1.108125  1.146921
## g_condition                            0.001806  0.005410  0.011624  0.029512
## g_customer                             0.001636  0.004531  0.009119  0.022339
## g_condition:customer                   0.002106  0.006912  0.015585  0.040756
##                                         97.5%
## mu                                    5.43998
## condition-control                     0.12432
## condition-experimental                0.01353
## customer-no                           0.09497
## customer-yes                          0.03476
## condition:customer-control.&.no       0.01152
## condition:customer-control.&.yes      0.12500
## condition:customer-experimental.&.no  0.12500
## condition:customer-experimental.&.yes 0.01152
## sig2                                  1.22542
## g_condition                           0.34480
## g_customer                            0.25701
## g_condition:customer                  0.52225

# Calculate posterior distributions for group means

# Main effects distributions
d.main.effects <- data.frame(chain.bayes.social.responsibility[, 2:5] + chain.bayes.social.responsibility[, 1])

# Interaction group distributions
d.control.yes <- chain.bayes.social.responsibility[, "mu"] + 
                 chain.bayes.social.responsibility[, "condition-control"] +
                 chain.bayes.social.responsibility[, "customer-yes"] + 
                 chain.bayes.social.responsibility[, "condition:customer-control.&.yes"]

d.control.no <- chain.bayes.social.responsibility[, "mu"] + 
                chain.bayes.social.responsibility[, "condition-control"] +
                chain.bayes.social.responsibility[, "customer-no"] + 
                chain.bayes.social.responsibility[, "condition:customer-control.&.no"]

d.experimental.yes <- chain.bayes.social.responsibility[, "mu"] + 
                      chain.bayes.social.responsibility[, "condition-experimental"] +
                      chain.bayes.social.responsibility[, "customer-yes"] + 
                      chain.bayes.social.responsibility[, "condition:customer-experimental.&.yes"]

d.experimental.no <- chain.bayes.social.responsibility[, "mu"] + 
                     chain.bayes.social.responsibility[, "condition-experimental"] +
                     chain.bayes.social.responsibility[, "customer-no"] + 
                     chain.bayes.social.responsibility[, "condition:customer-experimental.&.no"]

# Combine all posterior distributions
estimations <- bind_cols(as.data.frame(d.main.effects), 
                         condition.control_customer.yes = d.control.yes, 
                         condition.control_customer.no = d.control.no, 
                         condition.experimental_customer.yes = d.experimental.yes, 
                         condition.experimental_customer.no = d.experimental.no)

# Calculate 95% credible intervals
credible.intervals.soc.resp <- t(apply(estimations, 2, quantile, probs = c(0.025, 0.5, 0.975)))

# Save credible intervals in dataframe
credible.intervals.soc.resp.df <- data.frame(
  group = rownames(credible.intervals.soc.resp),
  median = credible.intervals.soc.resp[, "50%"],
  lower.ci = credible.intervals.soc.resp[, "2.5%"],
  upper.ci = credible.intervals.soc.resp[, "97.5%"]
)

rownames(credible.intervals.soc.resp.df) <- NULL
credible.intervals.soc.resp.df

The convergence diagnostics indicate that the model estimations were stable and converged appropriately (homogeneous patterns in the trace plots and normal distributions of the simulation results). The estimated means and credible intervals are consistent with the descriptive data. For better interpretation, we visualize them below.

# Visualization of posterior distributions
ggplot(credible.intervals.soc.resp.df [c(5:8), ], aes(x = group, y = median)) +
  
  # Point estimates with credible intervals
  geom_pointrange(aes(ymin = lower.ci, ymax = upper.ci), 
                  size = 0.7, color = "black") +
  
  # Point formatting for emphasis
  geom_point(size = 2, color = "black") +
  
  # Labels and titles
  labs(title = "Posterior Distributions - Social Responsibility Perception",
       subtitle = "Points show median estimates, bars show 95% credible intervals",
       y = "Social Responsibility Score (1-7 scale)", 
       x = "Experimental Group") +
  
  # Theme customization
  theme_bw() +
  coord_flip() +
  theme(plot.title = element_text(face = "bold"),
        plot.subtitle = element_text(color = "gray40"))

The posterior distributions provide evidence of the interaction we have been exploring: ControvPay differentially affects existing customers versus non-customers:

• Non-customers: Negligible difference (0.01 points) between those exposed and not exposed to ControvPay
• Existing customers: Decrease of 0.21 points in social responsibility perception when exposed to ControvPay

While the effect size is small, it could be meaningful given the large customer base. To further investigate this difference, we specifically evaluated the contrast between experimental conditions within the customer group using Bayesian post-hoc comparisons.

# Calculate posterior distribution of differences within customers
dif.in.customers <- d.control.yes - d.experimental.yes

# Calculate credible intervals for the difference
ci.dif.in.customers <- quantile(dif.in.customers, probs = c(0.025, 0.5, 0.975))
ci.dif.in.customers

##       2.5%        50%      97.5% 
## 0.02514849 0.20599950 0.41269383

The results confirm that social responsibility perception is significantly lower for the retailer’s customers exposed to ControvPay compared to those not exposed [median difference = 0.21 points, 95% CI: 0.03, 0.40].

The emergence of this difference specifically among existing customers aligns with theoretical expectations, as the retailer is known for high ethical standards, and its customers are likely more sensitive to ethical controversies in payment methods.

Given that social responsibility perception was associated with impactful outcomes like customer loyalty and purchase intention, we conducted mediation analyses within the customer subgroup to examine whether the effect on social responsibility had downstream consequences.

11. Path Analyses

Although prior analyses indicated no significant direct differences in loyalty and purchase intention between experimental conditions, it is possible that these effects were small and masked by statistical noise. Isolating the variability explained by indirect effects of the experimental condition on these variables through social responsibility perceptions can provide a more complete picture and more sensitivity for the experimental effects.

We used mediation analysis to examine these indirect effects with the R extension PROCESS (https://www.processmacro.org/index.html).

11.1. Checking the Fit of Path Analysis Models

We know collinearity between experimental conditions and social responsibility perception is not problematic for predicting other variables, as they are only slightly related (d < .30) (see Section 9.2).

To identify potential issues with outliers, normality, linearity, homoscedasticity, and independence of residuals, we:

• Specified the linear regression models implied in the path analysis.
• Plotted residual distributions and identified outliers’ influence on model estimates.

# Select subsample of retailer's customers
customers.pay.df <- pay.df[pay.df$customer == "yes", ]

# Model the two steps of path analysis for both DVs with multiple regression
m.path.1 <- lm(social.responsibility ~ condition, data = customers.pay.df)
m.path.2.loyalty <- lm(loyalty ~ condition + social.responsibility, data = customers.pay.df)
m.path.2.purchase.intention <- lm(purchase.intention ~ condition + social.responsibility, data = customers.pay.df)

# Obtain diagnostic plots
par(mfrow = c(2, 2))  # Set layout for 2×2 diagnostic plots
plot(m.path.1)

plot(m.path.2.loyalty) 

plot(m.path.2.purchase.intention)

In the residual vs. fitted values plots, we observe homogeneous patterns across fitted values, suggesting linear relationships and absence of heteroscedasticity issues.

The normal Q-Q plots indicate that standardized residuals closely follow the theoretical normal distribution, suggesting no major normality violations.

The residual vs. leverage plots identify potentially influential outliers based on high standardized residuals and leverage (cases 606, 786, 179, 294, 310, 23, 557, and 134). We will inspect these cases to identify potential data quality issues.

# Inspect potentially influential cases
customers.pay.df[customers.pay.df$id %in% c(606, 786, 179, 294, 310, 23, 557, 134), ]

We did not observe major data quality problems. Although some cases showed extreme scores in purchase intention or loyalty, these scores were consistent with other important variables like ease perceptions.

We therefore maintained all cases in the analyses to avoid artificial data manipulation.

11.2. Path Modeling for Purchase Intention

To analyze the mediation effect of experimental condition on purchase intention through social responsibility perceptions, we used the simple mediation model (Model 4) in PROCESS.

# Convert condition to numeric for PROCESS
customers.pay.df$condition_num <- ifelse(customers.pay.df$condition == "experimental", 1, 0)

# Run mediation analysis

## ATTENTION: Before running this function it is important to load PROCESS, 
## then remove "#" below for the function to run properly

#process(data = customers.pay.df, 
 #        y = "purchase.intention", 
 #        x = "condition_num", 
 #        m = "social.responsibility", 
 #        model = 4, 
 #        boot = 5000,
 #        seed = 123) 

To facilitate interpretation, we represent the results in a path diagram.

# Create path diagram with DOT notation, using the package DiagrammeR

grViz("
digraph mediation {
  graph [layout = neato, overlap = false, fontname = Arial]
  
  node [shape = rectangle, style = filled, fillcolor = 'lightblue', fontname = Arial]
  
  Condition [pos = '0,1!', label = 'Adding ControvPay']
  SocialResponsibility [pos = '2,2!', label = 'Social Responsibility\\nPerception']
  PurchaseIntention [pos = '4,1!', label = 'Purchase Intention']
  
  // Direct Effects
  Condition -> SocialResponsibility [label = 'a = -0.28*', fontsize = 12, color = 'black']
  SocialResponsibility -> PurchaseIntention [label = 'b = 0.09 (ns)', fontsize = 12, color = 'grey', style = 'dashed']
  Condition -> PurchaseIntention [label = 'c = 0.00 (ns)', fontsize = 12, color = 'grey', style = 'dashed']
  
  // Indirect Effect
  Indirect [shape = plaintext, pos = '2,0!', fillcolor = 'none']
  Indirect [label = 'Indirect Effect: a×b = -0.02\\n[95% CI: -0.07, 0.01]', fontsize = 14, color = 'darkred']
}
")

The results suggest that social responsibility perceptions do not translate into significant changes in immediate purchase intentions among existing customers.

This indirect effect is negligible (effect = -0.02) and its 95% confidence interval includes zero [-0.07, 0.01], indicating statistical non-significance.

11.3. Path Modeling for Loyalty

Following the same approach as the previous section, we used the simple mediation model (Model 4) in PROCESS to analyze the mediation effect of experimental condition on loyalty through social responsibility perceptions.

## ATTENTION: Before running this function it is important to load PROCESS, 
## then remove "#" below for the function to run properly

#process(data = customers.pay.df, 
 #        y = "loyalty", 
 #        x = "condition_num", 
 #        m = "social.responsibility", 
 #        model = 4, 
 #        boot = 5000,
 #        seed = 123) 

To facilitate visualization of the results, we represent them in a path diagram.

# Create path diagram with DOT notation
grViz("
digraph mediation {
  graph [layout = neato, overlap = false, fontname = Arial]
  
  node [shape = rectangle, style = filled, fillcolor = 'lightblue', fontname = Arial]
  
  Condition [pos = '0,1!', label = 'Adding ControvPay']
  SocialResponsibility [pos = '2,2!', label = 'Social Responsibility\\nPerception']
  Loyalty [pos = '4,1!', label = 'Loyalty']
  
  // Direct Effects
  Condition -> SocialResponsibility [label = 'a = -0.28*', fontsize = 12, color = 'black']
  SocialResponsibility -> Loyalty [label = 'b = 0.54***', fontsize = 12, color = 'black']
  Condition -> Loyalty [label = 'c = -0.03 (ns)', fontsize = 12, color = 'grey', style = 'dashed']
  
  // Indirect Effect
  Indirect [shape = plaintext, pos = '2,0!', fillcolor = 'none']
  Indirect [label = 'Indirect Effect: a×b = -0.15*\\n[95% CI: -0.27, -0.04]', fontsize = 14, color = 'darkred']
}
")

The path analysis confirms that social responsibility perceptions significantly mediate the relationship between ControvPay exposure and customer loyalty, accounting for a reduction of 0.15 points in the loyalty scale [95% CI: -0.27, -0.04]. We proceed with further analysis of this mediation effect’s practical implications in the following section.

12. Business Impact

Our analyses thus far indicate that introducing ControvPay may negatively affect social responsibility perceptions and loyalty among the retailer’s customers. To inform decision-making, we quantified the potential impact by estimating how many customers might transition from positive to negative perceptions or from loyal to non-loyal status.

For these predictions, we used the neutral points of the scales as cut-points:

• Social responsibility scores below 4 indicate negative perceptions
• Loyalty scores below 4 indicate non-loyal customers

To facilitate metric calculations, we created dichotomous numeric variables for both outcomes and the experimental condition.

# Create dichotomous variable for negative social responsibility perception
customers.pay.df$neg.social.responsibility <- NA_real_
customers.pay.df$neg.social.responsibility[customers.pay.df$social.responsibility < 4] <- 1
customers.pay.df$neg.social.responsibility[customers.pay.df$social.responsibility >= 4] <- 0

# Create dichotomous variable for non-loyal customers
customers.pay.df$unloyalty <- NA_real_
customers.pay.df$unloyalty[customers.pay.df$loyalty < 4] <- 1
customers.pay.df$unloyalty[customers.pay.df$loyalty >= 4] <- 0

# Create dichotomous variable for experimental condition
customers.pay.df$controvpay <- NA_real_
customers.pay.df$controvpay[customers.pay.df$condition == "experimental"] <- 1
customers.pay.df$controvpay[customers.pay.df$condition == "control"] <- 0

# Check distributions
summary(customers.pay.df)

##        id             condition        age          gender   
##  2      :  1   control     :192   Min.   :16.00   man  :197  
##  3      :  1   experimental:212   1st Qu.:32.00   woman:207  
##  8      :  1                      Median :41.00              
##  10     :  1                      Mean   :42.58              
##  13     :  1                      3rd Qu.:50.25              
##  14     :  1                      Max.   :85.00              
##  (Other):398                                                 
##           freq.customer customer     adoption      social.responsibility
##  none            :  0   no :  0   Min.   :0.0000   Min.   :2.500        
##  1 time          :162   yes:404   1st Qu.:0.0000   1st Qu.:4.500        
##  2-4 times       :110             Median :0.0000   Median :5.500        
##  5-9 times       : 89             Mean   :0.1179   Mean   :5.324        
##  10 times or more: 43             3rd Qu.:0.0000   3rd Qu.:6.000        
##                                   Max.   :1.0000   Max.   :7.000        
##                                   NA's   :192                           
##       ease      purchase.intention    loyalty              condition.customer
##  Min.   :1.00   Min.   :1.500      Min.   :1.000   ct.customers     :192     
##  1st Qu.:3.25   1st Qu.:4.000      1st Qu.:3.000   ct.non.customers :  0     
##  Median :4.25   Median :5.000      Median :4.000   exp.customers    :212     
##  Mean   :4.22   Mean   :4.984      Mean   :4.188   exp.non.customers:  0     
##  3rd Qu.:5.25   3rd Qu.:5.500      3rd Qu.:5.000                             
##  Max.   :7.00   Max.   :7.000      Max.   :7.000                             
##                                                                              
##  condition_num    neg.social.responsibility   unloyalty        controvpay    
##  Min.   :0.0000   Min.   :0.00000           Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.00000           1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :1.0000   Median :0.00000           Median :0.0000   Median :1.0000  
##  Mean   :0.5248   Mean   :0.07921           Mean   :0.3738   Mean   :0.5248  
##  3rd Qu.:1.0000   3rd Qu.:0.00000           3rd Qu.:1.0000   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :1.00000           Max.   :1.0000   Max.   :1.0000  
## 

The summary results indicate that 8% of the retailer’s customers reported negative social responsibility perceptions, and 37% reported disloyalty.

To predict how these proportions would vary with exposure to ControvPay, we used logistic regression models and a custom function from our toolkit to obtain interpretable predictive metrics. The function is specified below:

# Function to obtain key interpretation metrics for binomial logistic models

# Function name: f.metrics.binomial
#
# Input parameters:
#   my.model: fitted glm binomial model
#   my.matrix: design matrix for prediction scenarios
#
# Output: dataframe with comprehensive interpretation metrics including:
#   predicted probabilities, probability differences, odds ratios, and relative probabilities

f.metrics.binomial <- function(my.model, my.matrix){
  
  #--------1. Preliminary Checks: 
  
  # Making sure the model is glm binomial
  if (!inherits(my.model, "glm") || my.model$family$family != "binomial") {
    stop("The function requires a binomial model.")
  }
  
  
  #--------2. Odds Ratios
  
  # Obtaining the coefficients and confidence intervals
  coefs <- coef(my.model)
  ci <- confint(my.model)
  intercept <- coefs[1]
  
  # Exponentiation of the coefficients and their CIs to obtain odds ratios
  or <- exp(coefs)
  or_ci <- exp(ci)
  
  
  #--------3. Relative Probabilities (Risk Ratios)
  
  # Calculating the base probability from the intercept
  p0 <- or[1] / (1 + or[1])
  
  # Obtaining relative probabilities from risk ratios and base probability
  rp <- or / ((1 - p0) + (p0 * or))
  rp_ci_low <- or_ci[, 1] / ((1 - p0) + (p0 * or_ci[, 1]))
  rp_ci_high <- or_ci[, 2] / ((1 - p0) + (p0 * or_ci[, 2]))
  
  
  #--------4. Predicted Proportions
  
  # Extract log-odds and their sd for the individual presence of the factors (using the specified matrix)
  pred <- predict(my.model, newdata = my.matrix, type = "link", se.fit = TRUE)

  #Calculate confidence intervals in log-odds
  lower_logit <- pred$fit - 1.96 * pred$se.fit
  upper_logit <- pred$fit + 1.96 * pred$se.fit

  #Specify a function to transform log-odds to probabilities
  f.inv_logit <- function(x) {1 / (1 + exp(-x))}
  
  #Apply the function to the log-odds and the limits of their ci
  predicted_prob <- f.inv_logit(pred$fit) # Predicted probabilities
  lower_predicted_prob <- f.inv_logit(lower_logit) # Lower limit of CI for predicted probabilities
  upper_predicted_prob <- f.inv_logit(upper_logit) # Upper limit of CI for predicted probabilities

  
  #--------5. Differences in Predicted Probabilities via Bootstrap

  # Simulate a distribution of log-odds, for example with 10000 cases
  set.seed(123)
  simulated_log_odds <- as.data.frame(matrix(
  rnorm(10000 * length(pred$fit), mean = pred$fit, sd = pred$se.fit),
  ncol = length(pred$fit), 
  byrow = TRUE # the data is simulated line by line, so the values goes to the column of the corresponding factor
  ))

  # Transform logg-odds to expected probabilities (using the function specified above)
  simulated_probs <- f.inv_logit(simulated_log_odds)

  # Calculate differences from the first column, which correspond to the intercept probabilities
  diffs <- apply(simulated_probs, 2, function(col){
  col - simulated_probs[[1]]
  })

  # Calculate the means of the differences and 95% CI with percentiles
  prob_diff <- colMeans(diffs)
  lower_prob_diff  <- apply(diffs, 2, quantile, probs = 0.025)
  upper_prob_diff <- apply(diffs, 2, quantile, probs = 0.975)
  
  
  #-------- 6. Combine Results
  
  results <- data.frame(
    Term = names(rp),
    Prediction = as.numeric(predicted_prob)*100,
    Lower_Prediction = as.numeric(lower_predicted_prob)*100,
    Upper_Prediction = as.numeric(upper_predicted_prob)*100,
    Predicted_Change = as.numeric(prob_diff)*100,
    Lower_Predicted_Change = as.numeric(lower_prob_diff)*100,
    Upper_Predicted_Change = as.numeric(upper_prob_diff)*100,
    Relative_Probability = as.numeric(rp),
    Lower_Relative_Probability = as.numeric(rp_ci_low),
    Upper_Relative_Probability = as.numeric(rp_ci_high),
    Odds_Ratio = as.numeric(or),
    Lower_Odds_Ratio = as.numeric(or_ci[,1]),
    Upper_Odds_Ratio = as.numeric(or_ci[,2])
)
  
  # Remove meaningless metrics for intercept
  results[1, c("Predicted_Change",
             "Lower_Predicted_Change",
             "Upper_Predicted_Change",
             "Relative_Probability",
             "Lower_Relative_Probability",
             "Upper_Relative_Probability",
             "Odds_Ratio",
             "Lower_Odds_Ratio",
             "Upper_Odds_Ratio")] <- NA_real_

  
  # Return the results
  return(results)
}

After implementing these specifications, we calculated the expected impact of adding ControvPay on the retailer’s social responsibility perceptions and loyalty.

12.1. Impact on Social Responsibility Perception

To generate predictions, we specified a binomial logistic model predicting negative social responsibility perception from experimental condition within our customer dataframe, then obtained predictive metrics.

# Specify the model
m.neg.social.responsibility <- glm(neg.social.responsibility ~ controvpay, 
                                   data = customers.pay.df, 
                                   family = "binomial")
summary(m.neg.social.responsibility)

## 
## Call:
## glm(formula = neg.social.responsibility ~ controvpay, family = "binomial", 
##     data = customers.pay.df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.8006     0.3105  -9.019   <2e-16 ***
## controvpay    0.5929     0.3864   1.534    0.125    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 223.68  on 403  degrees of freedom
## Residual deviance: 221.22  on 402  degrees of freedom
## AIC: 225.22
## 
## Number of Fisher Scoring iterations: 5

# Specify the design matrix for predictions

# Extract factor names from the model
f.names.interc <- names(coef(m.neg.social.responsibility))
f.names <- f.names.interc[-1]

# Create design matrix with intercept and factor levels
my.matrix.df <- as.data.frame(diag(length(f.names.interc)))
my.matrix.df[, 1] <- f.names.interc
names(my.matrix.df) <- c("case", f.names)

# Check matrix structure
my.matrix.df

# Calculate model metrics
metrics.neg.social.responsibility <- f.metrics.binomial(my.model = m.neg.social.responsibility, 
                                                       my.matrix = my.matrix.df)

## Waiting for profiling to be done...

metrics.neg.social.responsibility

We observe that introducing ControvPay is associated with a 4% increase in customers holding negative perceptions of the retailer’s social responsibility within the customer base. This comparison is also displayed graphically in Section 12.3.

12.2. Impact on Loyalty

We followed a similar procedure to generate predictions for loyalty: * Specified a model predicting disloyalty from ControvPay exposure. * Created a design matrix representing cases with and without ControvPay exposure. * Obtained predictive metrics from the model and matrix.

# Specify the model
m.unloyalty <- glm(unloyalty ~ controvpay, 
                   data = customers.pay.df, 
                   family = "binomial")
summary(m.unloyalty)

## 
## Call:
## glm(formula = unloyalty ~ controvpay, family = "binomial", data = customers.pay.df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.6236     0.1514  -4.119 3.81e-05 ***
## controvpay    0.2024     0.2065   0.980    0.327    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 534.03  on 403  degrees of freedom
## Residual deviance: 533.07  on 402  degrees of freedom
## AIC: 537.07
## 
## Number of Fisher Scoring iterations: 4

# Specify the design matrix for predictions

# Extract factor names from the model
f.names.interc <- names(coef(m.unloyalty))
f.names <- f.names.interc[-1]

# Create design matrix
my.matrix.df <- as.data.frame(diag(length(f.names.interc)))
my.matrix.df[, 1] <- f.names.interc
names(my.matrix.df) <- c("case", f.names)

# Check matrix structure
my.matrix.df

# Calculate model metrics
metrics.unloyalty <- f.metrics.binomial(my.model = m.unloyalty, 
                                       my.matrix = my.matrix.df)

## Waiting for profiling to be done...

metrics.unloyalty

The introduction of ControvPay is expected to be associated with a 4.6% decrease in loyal customers within the retailer’s customer base. This comparison is also displayed graphically below.

12.3. Visual Representation of the Impact

We used the following code to visualize the impact of introducing ControvPay on the retailer’s customer base, examining both social responsibility perceptions and loyalty.

# Merge results
metrics.unloyalty$dep.variable <- "loyalty"
metrics.neg.social.responsibility$dep.variable <- "social.responsibility"
metrics.df <- bind_rows(metrics.neg.social.responsibility, metrics.unloyalty)

# Set DV order for logical presentation
metrics.df$dep.variable <- factor(metrics.df$dep.variable, 
                                 levels = c("social.responsibility", "loyalty"))

# Custom titles for faceted graphs
custom_titles <- c(
  social.responsibility = "Negative Perception of\nSocial Responsibility",
  loyalty = "Customer Disloyalty\n(Unlikely to Return)"
)

# Bar graph visualization
ggplot(metrics.df, aes(x = Term, y = Prediction)) +
  
  # Bar geometry
  geom_col(position = position_dodge(width = 0.8), width = 0.7, fill = "grey80") + 
  
  # Error bars for uncertainty
  geom_errorbar(
    aes(ymin = Lower_Prediction, ymax = Upper_Prediction),
    position = position_dodge(width = 0.8),
    width = 0.2) + 
  
  # Faceting by outcome variable
  facet_wrap(~ dep.variable, 
             labeller = labeller(dep.variable = custom_titles),
             strip.position = "top") + 
  
  # Labels and titles
  labs(
    y = "% of Retailer's Customers", 
    title = "Expected Impact of Adding ControvPay") + 
  
  # X-axis labels
  scale_x_discrete(labels = c("Current\nPayment Options", "With\nControvPay")) + 
  
  # Y-axis limits
  ylim(0, 60) + 
  
  # Data labels
  geom_text(
    aes(y = Upper_Prediction, label = paste0(round(Prediction, 0), "%")),
    position = position_dodge(width = 0.8),
    vjust = -1,
    size = 3.5
  ) +
  
  # Theme customization
  theme_minimal(base_size = 12) +
  theme(
    axis.title.x = element_blank(),
    strip.placement = "outside",
    panel.spacing = unit(1, "lines"),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    strip.text = element_text(face = "bold", size = 11),
    axis.text.x = element_text(angle = 0, vjust = 0, hjust = 0.5, size = 11),
    plot.margin = margin(10, 5, 10, 5)
  )

The graph illustrates the expected changes from introducing ControvPay among the retailer’s customers:

• Increased negative perceptions of social responsibility
• Increased customer disloyalty (decreased likelihood of website revisitation)

Although these predictions are subject to uncertainty, and the confidence intervals indicate the possibility of no change, the most probable outcome is the negative impact reflected by the point estimates.

13. Discussion

13.1. Managerial Implications

Our analysis suggests that although ControvPay is the preferred payment option for 3% of customers, its introduction is expected to yield primarily negative business consequences:

• It is not expected to increase purchase intention or perceived ease of the payment process
• Within the retailer’s customer base, it is projected to:
  • Increase negative perceptions of the retailer’s social responsibility by 4%
  • Increase customer disloyalty (decreased likelihood of return) by 4.6%

Given these risks and the absence of positive implications, we recommend deferring ControvPay implementation until further evidence demonstrates more favorable outcomes.

13.2. Limitations

Several limitations should be acknowledged. First, the sample size may have been insufficient to detect small effects, potentially limiting our power to identify other relevant impacts on purchase intention or additional outcomes where we found no statistical differences. Second, the analyses rely on self-reported data collected in an artificial environment—payment options presented in isolation. Evaluations using observational measures of real website interactions might yield different patterns of results.

13.3. Next Steps

Future research could extend these analyses by:

• Evaluating ControvPay’s impact through A/B testing of real customer interactions on the retailer’s website
• Conducting replication studies with larger sample sizes
• Exploring heterogeneity across customer segments to identify groups more sensitive to ControvPay introduction

14. Bibliography

• Maignan, I. (2001). Consumers’ perceptions of corporate social responsibilities: A cross-cultural comparison. Journal of business ethics, 30(1), 57-72.

• Sauro, J. (2015). SUPR-Q: A comprehensive measure of the quality of the website user experience. Journal of usability studies, 10(2).

• Venkatesh, V., & Davis, F. D. (2000). A Theoretical Extension of the Technology Acceptance Model: Four Longitudinal Field Studies. Management Science, 46(2), 186-204. https://doi.org/10.1287/mnsc.46.2.186.11926

Appendix

To save the dataframes we can use the code below

write.csv(pay.df, "pay.df.csv")