概要
このノートは、高松市ウェルビーイング調査データを使った
lvsem のデモです。
wellbeing.csv
は、1行目に領域名、2行目に質問文、3行目以降に回答データが入っています。
ここでは質問項目を5つのブロックにまとめ、Wellbeing
を結果変数とする構造モデルを推定します。
分析の流れは次の通りです。
- CSVを読み込み、項目と領域を整理する
- 5つの分析ブロックを定義する
lvsem で潜在変数スコアを作成する- 構造パスと効果を推定する
- 測定モデルの信頼性、AVE、外部負荷量を確認する
- モデル図を作成する
パッケージ
rm(list = ls(all = TRUE))
library(lvsem)
library(ggplot2)
library(knitr)
データ読み込み
data_file <- "wellbeing.csv"
data_candidates <- c(
file.path("data000", data_file),
data_file,
file.path("..", "data000", 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)
meta <- read.csv(
data_path,
header = FALSE,
nrows = 2,
check.names = FALSE,
fileEncoding = "CP932"
)
category_jp <- as.character(meta[1, -1])
item_label <- as.character(meta[2, -1])
raw_data <- read.csv(
data_path,
header = FALSE,
skip = 2,
check.names = FALSE,
fileEncoding = "CP932"
)
dim(raw_data)
## [1] 100 56
head(raw_data)
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20 V21
## 1 10001 NA 6 NA 5 6 4 4 4 4 4 5 2 3 4 3 1 2 2
## 2 10002 NA 7 NA 7 7 4 4 4 4 4 4 2 4 4 4 4 4 4
## 3 10003 NA 10 NA 9 9 4 5 5 5 5 5 1 2 4 4 4 5 5
## 4 10004 NA 5 NA 5 5 3 5 3 3 1 3 1 4 1 2 4 4 4
## 5 10005 NA 5 NA 7 6 3 4 3 4 3 3 2 4 3 2 4 3 4
## 6 10006 NA 7 NA 9 6 4 4 4 4 4 5 1 4 3 3 3 4 4
## V22 V23 V24 V25 V26 V27 V28 V29 V30 V31 V32 V33 V34 V35 V36 V37 V38 V39 V40
## 1 4 4 4 2 4 4 4 4 4 4 5 2 4 4 4 5 3 4 4
## 2 4 4 4 4 4 4 3 4 4 4 4 4 4 4 4 4 4 4 4
## 3 5 4 4 4 4 5 4 4 3 4 4 4 4 5 5 4 3 3 3
## 4 4 4 4 3 3 4 3 2 3 4 2 3 4 2 2 3 3 3 4
## 5 3 4 3 3 3 4 3 4 4 3 3 3 3 4 4 4 3 3 4
## 6 4 4 4 3 4 4 4 4 4 4 4 3 3 3 4 4 3 3 4
## V41 V42 V43 V44 V45 V46 V47 V48 V49 V50 V51 V52 V53 V54 V55 V56
## 1 4 4 2 3 4 3 3 2 2 3 3 2 2 4 3
## 2 4 4 4 4 4 3 4 4 4 4 4 4 4 4 4
## 3 4 3 2 3 4 4 4 4 5 5 5 5 4 5 4
## 4 3 3 2 1 3 3 4 3 4 3 4 4 3 3 3
## 5 4 3 2 4 3 3 3 3 2 3 4 3 3 3 3
## 6 4 3 3 3 4 4 4 4 5 3 4 4 4 4 4
領域の整理
category_keys <- c(
"attr", "wellBeing", "healthcare", "shoppingDining", "housing",
"mobility", "recreation", "childcare", "primarySecondaryEdu",
"localGovernment", "digitalLife", "publicSpace", "urbanLandscape",
"naturalLandscape", "naturalResources", "environmentalHarmony",
"naturalDisasters", "accidentsCrime", "communityConnectedness",
"diversityTolerance", "selfEfficacy", "healthStatus", "cultureArts",
"educationalOpportunities", "employmentIncome", "entrepreneurship",
"original"
)
category_order <- unique(category_jp)
category_index <- match(category_jp, category_order)
item_keys <- make.unique(category_keys[category_index], sep = "_")
names(raw_data) <- c("RespondentID", item_keys)
item_dictionary <- data.frame(
category_key = category_keys[category_index],
category = category_jp,
item_key = item_keys,
item = item_label,
stringsAsFactors = FALSE
)
kable(head(item_dictionary, 12), caption = "項目辞書の先頭")
項目辞書の先頭
| attr |
属性 |
attr |
性別 |
| attr |
属性 |
attr_1 |
年代 |
| attr |
属性 |
attr_2 |
住所 |
| wellBeing |
幸福度・満足度 |
wellBeing |
現在、あなたはどの程度幸せですか。「とても幸せ」を 10
点、「とても不幸」を0点とすると、何点くらいになると思いますか。いずれかの数字を1つだけお答えください。 |
| wellBeing |
幸福度・満足度 |
wellBeing_1 |
今から5年後、あなたはどの程度幸せだと思いますか。「とても幸せ」を
10
点、「とても不幸」を0点とすると、何点くらいになると思いますか。いずれかの数字を1つだけお答えください。 |
| wellBeing |
幸福度・満足度 |
wellBeing_2 |
現在、あなたの住んでいる地域の暮らしにどの程度満足していますか。「とても満足」を
10
点、「とても不満足」を0点とすると、何点くらいになると思いますか。いずれかの数字を1つだけお答えください。 |
| wellBeing |
幸福度・満足度 |
wellBeing_3 |
あなたの町内（集落）の人々は、大体において、どれぐらい幸せだと思いますか。「とても幸せ」を10点、「とても不幸」を0点として、いずれかの数字を1つだけお答えください。ここでは自分の同居家族は除いて考えてください。 |
| wellBeing |
幸福度・満足度 |
wellBeing_4 |
自分だけでなく、身近なまわりの人も楽しい気持ちでいると思う |
| healthcare |
医療・福祉 |
healthcare |
暮らしている地域は、医療機関が充実している |
| healthcare |
医療・福祉 |
healthcare_1 |
私の暮らしている地域では、介護・福祉施設のサービスが受けやすい |
| shoppingDining |
買物・飲食 |
shoppingDining |
暮らしている地域は、日常の買い物にまったく不便がない |
| shoppingDining |
買物・飲食 |
shoppingDining_1 |
私の暮らしている地域では、飲食を楽しめる場所が充実している |
モデル定義
このデモでは、調査項目を次の5ブロックにまとめます。
Wellbeing: 幸福度・満足度Health: 健康状態、医療・福祉Social: 地域とのつながりEnv:
住宅環境、移動・交通、公共空間、自然景観、都市景観EcoCult:
雇用・所得、文化・芸術、教育機会の豊かさ、買物・飲食
grouped_items <- split(item_keys, category_keys[category_index])
blocks <- list(
Wellbeing = grouped_items[["wellBeing"]],
Health = c(grouped_items[["healthStatus"]], grouped_items[["healthcare"]]),
Social = grouped_items[["communityConnectedness"]],
Env = c(
grouped_items[["housing"]],
grouped_items[["mobility"]],
grouped_items[["publicSpace"]],
grouped_items[["naturalLandscape"]],
grouped_items[["urbanLandscape"]]
),
EcoCult = c(
grouped_items[["employmentIncome"]],
grouped_items[["cultureArts"]],
grouped_items[["educationalOpportunities"]],
grouped_items[["shoppingDining"]]
)
)
block_summary <- data.frame(
block = names(blocks),
item_count = vapply(blocks, length, integer(1))
)
kable(block_summary, caption = "分析ブロック")
分析ブロック
| Wellbeing |
Wellbeing |
5 |
| Health |
Health |
4 |
| Social |
Social |
5 |
| Env |
Env |
8 |
| EcoCult |
EcoCult |
7 |
latent_vars <- names(blocks)
path <- matrix(
0,
nrow = length(latent_vars),
ncol = length(latent_vars),
dimnames = list(latent_vars, latent_vars)
)
path["Wellbeing", c("Health", "Social", "Env", "EcoCult")] <- 1
path
## Wellbeing Health Social Env EcoCult
## Wellbeing 0 1 1 1 1
## Health 0 0 0 0 0
## Social 0 0 0 0 0
## Env 0 0 0 0 0
## EcoCult 0 0 0 0 0
score_relation <- path + t(path)
score_relation
## Wellbeing Health Social Env EcoCult
## Wellbeing 0 1 1 1 1
## Health 1 0 0 0 0
## Social 1 0 0 0 0
## Env 1 0 0 0 0
## EcoCult 1 0 0 0 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
})
too_constant <- vapply(model_data, function(x) {
observed <- x[!is.na(x)]
if (!length(observed)) {
return(TRUE)
}
stats::sd(observed) == 0 || max(table(observed)) / length(observed) >= 0.80
}, logical(1))
if (any(too_constant)) {
model_data <- model_data[, !too_constant, drop = FALSE]
blocks <- lapply(blocks, function(x) intersect(x, names(model_data)))
}
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),
removed_constant_variables = sum(too_constant)
)
## samples variables removed_constant_variables
## 1 100 28 1
潜在変数スコア
lv_scores <- make_lv_scores(model_data, blocks, score_relation)
head(lv_scores)
## Wellbeing Health Social Env EcoCult
## 1 -1.1554429 -0.90189036 0.45910930 0.6250015 -0.6699975
## 2 0.1611687 0.23998621 0.51821763 0.1330351 0.9992083
## 3 2.3606732 0.95508628 0.07456707 0.9607024 2.5197512
## 4 -2.3980084 0.05001247 -0.19069374 -1.0766676 -0.6876581
## 5 -1.3601242 -0.77212340 0.13997187 -0.3311867 -0.2668526
## 6 0.5481766 0.51408977 0.13997187 0.5891602 0.7410003
round(cor(lv_scores), 3)
## Wellbeing Health Social Env EcoCult
## Wellbeing 1.000 0.674 0.440 0.716 0.704
## Health 0.674 1.000 0.211 0.724 0.654
## Social 0.440 0.211 1.000 0.324 0.323
## Env 0.716 0.724 0.324 1.000 0.684
## EcoCult 0.704 0.654 0.323 0.684 1.000
構造モデル
path_coefficients <- estimate_paths(lv_scores, path, scale_scores = FALSE)
path_effects <- effects_from_beta(path_coefficients)
path_coefficients
## path beta
## 1 Health_to_Wellbeing 0.6678333
## 2 Social_to_Wellbeing 0.6564387
## 3 Env_to_Wellbeing 0.4868310
## 4 EcoCult_to_Wellbeing 0.4247008
path_effects
## relationships direct indirect total
## 1 EcoCult -> Wellbeing 0.4247008 0 0.4247008
## 2 Env -> Wellbeing 0.4868310 0 0.4868310
## 3 Health -> Wellbeing 0.6678333 0 0.6678333
## 4 Social -> Wellbeing 0.6564387 0 0.6564387
lm_fit <- lm(
Wellbeing ~ Health + Social + Env + EcoCult,
data = lv_scores
)
structural_table <- data.frame(
variable = rownames(summary(lm_fit)$coefficients),
summary(lm_fit)$coefficients,
row.names = NULL,
check.names = FALSE
)
kable(
data.frame(
variable = structural_table$variable[-1],
round(structural_table[-1, -1], 3)
),
caption = "構造モデルの推定結果"
)
構造モデルの推定結果
| 2 |
Health |
0.668 |
0.258 |
2.592 |
0.011 |
| 3 |
Social |
0.656 |
0.206 |
3.183 |
0.002 |
| 4 |
Env |
0.487 |
0.170 |
2.864 |
0.005 |
| 5 |
EcoCult |
0.425 |
0.129 |
3.291 |
0.001 |
r2_wellbeing <- summary(lm_fit)$r.squared
adj_r2_wellbeing <- summary(lm_fit)$adj.r.squared
r2_wellbeing
## [1] 0.6550963
adj_r2_wellbeing
## [1] 0.6405741
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による構造パス
| Health_to_Wellbeing |
0.239 |
0.110 |
2.171 |
| Social_to_Wellbeing |
0.206 |
0.088 |
2.329 |
| Env_to_Wellbeing |
0.277 |
0.119 |
2.331 |
| EcoCult_to_Wellbeing |
0.291 |
0.102 |
2.855 |
測定モデルの要約
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
| Wellbeing |
0.808 |
0.876 |
0.641 |
| Health |
0.530 |
0.740 |
0.433 |
| Social |
0.761 |
0.830 |
0.496 |
| Env |
0.640 |
0.749 |
0.339 |
| EcoCult |
0.804 |
0.856 |
0.462 |
kable(
data.frame(
block = measures$loadings$block,
item = measures$loadings$item,
loading = round(measures$loadings$loading, 3)
),
caption = "外部負荷量"
)
外部負荷量
| Wellbeing |
wellBeing |
0.791 |
| Wellbeing |
wellBeing_2 |
0.888 |
| Wellbeing |
wellBeing_3 |
0.847 |
| Wellbeing |
wellBeing_4 |
0.659 |
| Health |
healthStatus |
0.713 |
| Health |
healthStatus_1 |
0.810 |
| Health |
healthcare |
0.668 |
| Health |
healthcare_1 |
0.349 |
| Social |
communityConnectedness |
0.595 |
| Social |
communityConnectedness_1 |
0.752 |
| Social |
communityConnectedness_2 |
0.777 |
| Social |
communityConnectedness_3 |
0.632 |
| Social |
communityConnectedness_4 |
0.748 |
| Env |
housing |
0.621 |
| Env |
housing_1 |
-0.194 |
| Env |
housing_2 |
0.219 |
| Env |
mobility |
0.537 |
| Env |
publicSpace |
0.720 |
| Env |
publicSpace_1 |
0.707 |
| Env |
naturalLandscape |
0.625 |
| Env |
urbanLandscape |
0.735 |
| EcoCult |
employmentIncome |
0.759 |
| EcoCult |
employmentIncome_1 |
0.655 |
| EcoCult |
cultureArts |
0.742 |
| EcoCult |
cultureArts_1 |
0.580 |
| EcoCult |
educationalOpportunities |
0.721 |
| EcoCult |
shoppingDining |
0.595 |
| EcoCult |
shoppingDining_1 |
0.686 |
モデル図
node_df <- data.frame(
block = c("Health", "Social", "Env", "EcoCult", "Wellbeing"),
label = c(
"健康・医療\n(Health)",
"地域とのつながり\n(Social)",
"住環境・自然\n(Env)",
"経済・文化\n(EcoCult)",
"幸福度・満足度\n(Wellbeing)"
),
x = c(0, 0, 0, 0, 3.2),
y = c(3, 2, 1, 0, 1.5)
)
edge_df <- data.frame(
from = sub("_to_Wellbeing", "", path_coefficients$path),
to = "Wellbeing",
beta = path_coefficients$beta,
label = sprintf("beta = %.3f", path_coefficients$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 + c(-0.18, 0.16, 0.10, -0.12)
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 <- ifelse(item_df$block == "Wellbeing", item_df$x_node + 1.25, item_df$x_node - 1.20)
item_df$y <- item_df$y_node + (item_df$i - (item_df$n + 1) / 2) * 0.15
item_df$label <- sprintf("%s\n%.3f", item_df$item, item_df$loading)
model_plot <- ggplot() +
geom_segment(
data = edge_df,
aes(x = x_from + 0.45, y = y_from, xend = x_to - 0.75, yend = y_to),
arrow = arrow(length = unit(0.22, "cm")),
linewidth = 0.8,
color = "#2f3a45"
) +
geom_label(
data = item_df,
aes(x = x, y = y, label = label),
size = 2.6,
label.size = 0.25,
fill = "#f7fbff",
color = "#2a3340"
) +
geom_segment(
data = item_df,
aes(
x = ifelse(block == "Wellbeing", x - 0.30, x + 0.30),
y = y,
xend = ifelse(block == "Wellbeing", x_node + 0.82, x_node - 0.62),
yend = y_node
),
linewidth = 0.35,
color = "#8b97a5"
) +
geom_label(
data = node_df,
aes(x = x, y = y, label = label),
size = 4,
fontface = "bold",
label.size = 0.35,
fill = "#ffffff",
color = "#111827"
) +
geom_label(
data = edge_df,
aes(x = x_label, y = y_label, label = label),
size = 3,
label.size = 0,
fill = "#ffffff",
color = "#111827"
) +
annotate(
"label",
x = 3.2,
y = 0.72,
label = sprintf("R2 = %.3f", r2_wellbeing),
size = 3.6,
label.size = 0.25,
fill = "#eef7f2"
) +
coord_cartesian(xlim = c(-1.75, 5.10), ylim = c(-0.55, 3.55), clip = "off") +
theme_void() +
labs(title = "高松市ウェルビーイングデータのlvsemモデル")
model_plot
ggsave(file.path(data_dir, "data000_lvsem_wellbeing_model.png"),
model_plot, width = 11, height = 7, dpi = 180)
出力の確認
cat("サンプルサイズ:", nrow(model_data), "\n")
## サンプルサイズ: 100
cat("ウェルビーイングのR2:", round(r2_wellbeing, 3), "\n")
## ウェルビーイングのR2: 0.655
cat("調整済みR2:", round(adj_r2_wellbeing, 3), "\n")
## 調整済みR2: 0.641
参考
- ローカルパッケージ:
lvsem
- データ:
wellbeing.csv