Data002: 乳製品ブランドデータ lvsem デモ

概要

このノートは、乳製品ブランドに関する公開データを使った lvsem のデモです。

データには、消費者満足、ブランド態度、ブランド選好、購買意図に関する回答が含まれています。 ここでは、論文で示されている構造に沿って、CS から BA、BP、PI へ、 さらに BA から BP と PI へ、BP から PI へ向かうモデルを推定します。

分析の流れは次の通りです。

• CSVを読み込む
• 観測変数ブロックと構造パス行列を定義する
• lvsem で潜在変数スコアを作成する
• パス係数、効果、R2を推定する
• bootstrapでパス係数のばらつきを確認する
• 測定モデルの信頼性、AVE、外部負荷量を確認する
• モデル図を作成する

パッケージ

rm(list = ls(all = TRUE))

library(lvsem)
library(ggplot2)
library(knitr)

データ読み込み

data_file <- "Data Base.csv"
data_candidates <- c(
  file.path("data002", data_file),
  data_file,
  file.path("..", "data002", data_file)
)
data_path <- data_candidates[file.exists(data_candidates)][1]

if (is.na(data_path)) {
  stop("CSVファイルが見つかりません: ", data_file)
}

data_dir <- dirname(data_path)

raw_data <- read.csv(data_path, sep = ";", check.names = TRUE)

dim(raw_data)

## [1] 195  19

head(raw_data)

##   Gender Age MaritalStatus CS1 CS2 CS3 CS4 CS5 BP1 BP2 BP3 BP4 BA1 BA2 BA3 PI1
## 1  women  19        single   5   5   5   5   5   5   3   3   5   5   5   5   5
## 2  women  19        single   4   5   5   5   4   5   5   5   5   5   5   5   5
## 3    men  24        single   4   4   4   4   3   1   2   2   1   4   5   3   3
## 4    men  44      married    5   5   5   5   5   3   3   2   2   4   4   4   4
## 5  women  21        single   3   4   3   4   3   1   1   1   1   3   4   5   3
## 6    men  31        single   5   5   5   5   4   4   5   5   5   5   5   5   5
##   PI2 PI3 PI4
## 1   5   5   5
## 2   5   5   5
## 3   4   3   4
## 4   3   3   3
## 5   2   3   3
## 6   4   4   4

モデル定義

このデモでは、次の4つの潜在変数を使います。

• CS: 消費者満足
• BA: ブランド態度
• BP: ブランド選好
• PI: 購買意図

blocks <- list(
  CS = paste0("CS", 1:5),
  BA = paste0("BA", 1:3),
  BP = paste0("BP", 1:4),
  PI = paste0("PI", 1:4)
)

latent_vars <- names(blocks)
path <- matrix(
  0,
  nrow = length(latent_vars),
  ncol = length(latent_vars),
  dimnames = list(latent_vars, latent_vars)
)

path["BA", "CS"] <- 1
path["BP", c("CS", "BA")] <- 1
path["PI", c("CS", "BA", "BP")] <- 1
path

##    CS BA BP PI
## CS  0  0  0  0
## BA  1  0  0  0
## BP  1  1  0  0
## PI  1  1  1  0

score_relation <- path + t(path)
score_relation

##    CS BA BP PI
## CS  0  1  1  1
## BA  1  0  1  1
## BP  1  1  0  1
## PI  1  1  1  0

前処理

model_vars <- unique(unlist(blocks))
model_data <- raw_data[, model_vars, drop = FALSE]

model_data[] <- lapply(model_data, function(x) {
  x <- as.numeric(as.character(x))
  x[!is.finite(x)] <- NA_real_
  x
})

model_data[] <- lapply(model_data, function(x) {
  ifelse(is.na(x), stats::median(x, na.rm = TRUE), x)
})

data.frame(
  samples = nrow(model_data),
  variables = ncol(model_data)
)

##   samples variables
## 1     195        16

潜在変数スコア

lv_scores <- make_lv_scores(model_data, blocks, score_relation)

head(lv_scores)

##            CS         BA         BP         PI
## 1  1.03805651  0.6105095  0.4050641  0.9846493
## 2  0.70608214  0.6105095  1.1150581  0.9846493
## 3 -0.06122641 -0.1121397 -1.3732591 -0.3688359
## 4  1.03805651 -0.1241319 -0.6619449 -0.5598812
## 5 -0.45650242 -0.1492832 -1.7282561 -1.0333737
## 6  0.89603355  0.6105095  0.9321093  0.3368683

round(cor(lv_scores), 3)

##       CS    BA    BP    PI
## CS 1.000 0.583 0.671 0.779
## BA 0.583 1.000 0.529 0.691
## BP 0.671 0.529 1.000 0.717
## PI 0.779 0.691 0.717 1.000

構造モデル

path_coefficients <- estimate_paths(lv_scores, path)
path_effects <- effects_from_beta(path_coefficients)

path_coefficients

##       path      beta
## 1 CS_to_BA 0.5827942
## 2 CS_to_BP 0.5492620
## 3 BA_to_BP 0.2089492
## 4 CS_to_PI 0.4165730
## 5 BA_to_PI 0.3002360
## 6 BP_to_PI 0.2788174

path_effects

##   relationships    direct   indirect     total
## 1      BA -> BP 0.2089492 0.00000000 0.2089492
## 2      BA -> PI 0.3002360 0.05825867 0.3584947
## 3      BP -> PI 0.2788174 0.00000000 0.2788174
## 4      CS -> BA 0.5827942 0.00000000 0.5827942
## 5      CS -> BP 0.5492620 0.12177438 0.6710363
## 6      CS -> PI 0.4165730 0.36207242 0.7786455

r2_values <- sapply(c("BA", "BP", "PI"), function(dv) {
  preds <- names(which(path[dv, ] == 1))
  fit <- lm(
    as.formula(paste(dv, "~", paste(preds, collapse = " + "))),
    data = as.data.frame(scale(lv_scores))
  )
  summary(fit)$r.squared
})

kable(
  data.frame(
    construct = names(r2_values),
    R2 = round(as.numeric(r2_values), 3)
  ),
  caption = "内生潜在変数のR2"
)

内生潜在変数のR2

construct	R2
BA	0.340
BP	0.479
PI	0.732

Bootstrap

set.seed(123)
boot_paths <- bootstrap_paths(
  model_data,
  blocks,
  path,
  br = 1000,
  seed = 123,
  score_matrix = score_relation
)

kable(
  data.frame(
    path = boot_paths$path,
    round(boot_paths[, c("beta", "boot_se", "t_value")], 3)
  ),
  caption = "Bootstrapによる構造パス"
)

Bootstrapによる構造パス

path	beta	boot_se	t_value
CS_to_BA	0.583	0.061	9.512
CS_to_BP	0.549	0.070	7.881
BA_to_BP	0.209	0.067	3.137
CS_to_PI	0.417	0.066	6.318
BA_to_PI	0.300	0.053	5.666
BP_to_PI	0.279	0.066	4.252

測定モデルの要約

measures <- measurement_summary(model_data, blocks, lv_scores)

kable(
  data.frame(
    block = measures$reliability$block,
    round(measures$reliability[, c("alpha", "composite_reliability", "ave")], 3)
  ),
  caption = "信頼性とAVE"
)

信頼性とAVE

block	alpha	composite_reliability	ave
CS	0.945	0.959	0.825
BA	0.927	0.953	0.872
BP	0.958	0.970	0.889
PI	0.886	0.923	0.751

kable(
  data.frame(
    block = measures$loadings$block,
    item = measures$loadings$item,
    loading = round(measures$loadings$loading, 3)
  ),
  caption = "外部負荷量"
)

外部負荷量

block	item	loading
CS	CS1	0.928
CS	CS2	0.939
CS	CS3	0.947
CS	CS4	0.942
CS	CS5	0.772
BA	BA1	0.935
BA	BA2	0.937
BA	BA3	0.929
BP	BP1	0.936
BP	BP2	0.943
BP	BP3	0.949
BP	BP4	0.943
PI	PI1	0.721
PI	PI2	0.925
PI	PI3	0.896
PI	PI4	0.910

モデル図

node_df <- data.frame(
  block = c("CS", "BA", "BP", "PI"),
  label = c(
    "消費者満足\n(CS)",
    "ブランド態度\n(BA)",
    "ブランド選好\n(BP)",
    "購買意図\n(PI)"
  ),
  x = c(0, 2.2, 2.2, 4.4),
  y = c(1.4, 2.2, 0.6, 1.4)
)

edge_df <- path_coefficients
edge_split <- strsplit(edge_df$path, "_to_")
edge_df$from <- vapply(edge_split, `[`, character(1), 1)
edge_df$to <- vapply(edge_split, `[`, character(1), 2)
edge_df$label <- sprintf("beta = %.3f", edge_df$beta)

edge_df <- merge(
  edge_df,
  node_df[, c("block", "x", "y")],
  by.x = "from",
  by.y = "block"
)
names(edge_df)[names(edge_df) %in% c("x", "y")] <- c("x_from", "y_from")

edge_df <- merge(
  edge_df,
  node_df[, c("block", "x", "y")],
  by.x = "to",
  by.y = "block"
)
names(edge_df)[names(edge_df) %in% c("x", "y")] <- c("x_to", "y_to")

edge_df$x_label <- (edge_df$x_from + edge_df$x_to) / 2
edge_df$y_label <- (edge_df$y_from + edge_df$y_to) / 2

item_df <- measures$loadings
item_df$i <- ave(seq_len(nrow(item_df)), item_df$block, FUN = seq_along)
item_df$n <- ave(seq_len(nrow(item_df)), item_df$block, FUN = length)
item_df <- merge(item_df, node_df[, c("block", "x", "y")], by = "block")
names(item_df)[names(item_df) %in% c("x", "y")] <- c("x_node", "y_node")
item_df$x <- item_df$x_node
item_df$y <- item_df$y_node - 0.70 - (item_df$i - 1) * 0.17
item_df$label <- sprintf("%s %.3f", item_df$item, item_df$loading)

r2_df <- data.frame(
  block = names(r2_values),
  r2_label = sprintf("R2 = %.3f", as.numeric(r2_values))
)
r2_df <- merge(r2_df, node_df[, c("block", "x", "y")], by = "block")

model_plot <- ggplot() +
  geom_segment(
    data = edge_df,
    aes(x = x_from + 0.45, y = y_from, xend = x_to - 0.45, yend = y_to),
    arrow = arrow(length = unit(0.22, "cm")),
    linewidth = 0.8,
    color = "#2f3a45"
  ) +
  geom_label(
    data = node_df,
    aes(x = x, y = y, label = label),
    size = 4.1,
    fontface = "bold",
    label.size = 0.35,
    fill = "#ffffff",
    color = "#111827"
  ) +
  geom_label(
    data = edge_df,
    aes(x = x_label, y = y_label + 0.16, label = label),
    size = 3,
    label.size = 0,
    fill = "#ffffff",
    color = "#111827"
  ) +
  geom_label(
    data = item_df,
    aes(x = x, y = y, label = label),
    size = 2.8,
    label.size = 0.2,
    fill = "#f7fbff",
    color = "#2a3340"
  ) +
  geom_label(
    data = r2_df,
    aes(x = x, y = y + 0.55, label = r2_label),
    size = 3.1,
    label.size = 0.25,
    fill = "#eef7f2",
    color = "#111827"
  ) +
  coord_cartesian(xlim = c(-0.9, 5.3), ylim = c(-0.7, 3.0), clip = "off") +
  theme_void() +
  labs(title = "乳製品ブランドデータのlvsemモデル")

model_plot

ggsave(file.path(data_dir, "data002_lvsem_brand_model.png"),
       model_plot, width = 11, height = 7, dpi = 180)

出力の確認

cat("サンプルサイズ:", nrow(model_data), "\n")

## サンプルサイズ: 195

print(round(r2_values, 3))

##    BA    BP    PI 
## 0.340 0.479 0.732

参考

• Boubker & Douayri (2020), Data in Brief, DOI: 10.1016/j.dib.2020.106172
• ローカルパッケージ: lvsem
• データ: Data Base.csv