本サイトは、ウィラワン ドニ ダハナ, 勝又 壮太郎 (2023).『Rによるマーケティング・データ分析』新世社のサポートサイトです。
更新履歴
データをまとめてダウンロードしたい方はこちら。
本文で解説しているコードは以下にあります。開いて確認してください。
#-2.1--------------------------------------------------------------------------#
#データのインポート
data_chap2 <- read.csv("chapter_02.csv", header = TRUE)
head(data_chap2)
#-2.2--------------------------------------------------------------------------#
# 店舗1の売上高
store1_sales <- data_chap2[,"Store1"]
# 算術平均
mean(store1_sales)
## [1] 37.96318
# 幾何平均
exp(mean(log(store1_sales)))
## [1] 37.38505
# 調和平均
length(store1_sales) / sum(1 / store1_sales)
## [1] 36.73915
#-2.3--------------------------------------------------------------------------#
# 分散
var(store1_sales)
## [1] 40.35984
# 標本分散、不偏分散を定義通りに計算する
N <- length(store1_sales)
m <- mean(store1_sales)
(1 / N) * sum((store1_sales - m)^2)
## [1] 40.24927
(1 / (N - 1)) * sum((store1_sales - m)^2)
## [1] 40.35984
#-2.4--------------------------------------------------------------------------#
hist(store1_sales)
#-2.5--------------------------------------------------------------------------#
summary(store1_sales)
#-2.6--------------------------------------------------------------------------#
hist_out <- hist(store1_sales)
hist_out$mids[which.max(hist_out$counts)]
#-2.7--------------------------------------------------------------------------#
weather <- data_chap2[,"Weather"]
table(weather)
#-2.8--------------------------------------------------------------------------#
store1_rank <- data_chap2[,"Store1Rank"]
table(store1_rank)
mean(store1_rank)
store2_rank <- data_chap2[,"Store2Rank"]
table(store2_rank)
mean(store2_rank)
#-2.9--------------------------------------------------------------------------#
table(store1_rank, weather)
table(store2_rank, weather)
#-2.10-------------------------------------------------------------------------#
chisq.test(table(store1_rank, weather))
chisq.test(table(store2_rank, weather))
#-2.11-------------------------------------------------------------------------#
store1_sunny <- store1_sales[weather == 1]
store1_rain <- store1_sales[weather == 2]
var(store1_sunny)
var(store1_rain)
var.test(store1_rain,store1_sunny)
#-2.12-------------------------------------------------------------------------#
t.test(store1_sunny, store1_rain, var.equal = TRUE)
#-2.13-------------------------------------------------------------------------#
rain <- data_chap2[,"Rain"]
plot(store1_sales, rain)
#-2.14-------------------------------------------------------------------------#
cor(store1_sales, rain)
#-2.15-------------------------------------------------------------------------#
mean(store1_sales)
#店舗2の売上オブジェクトstore2_salesを定義します。
store2_sales <- data_chap2[,"Store2"]
mean(store2_sales)
#店舗3の売上オブジェクトstore3_salesを定義します。
store3_sales <- data_chap2[,"Store3"]
mean(store3_sales)
#-2.16-------------------------------------------------------------------------#
all_sales <- c(store1_sales, store2_sales, store3_sales)
#-2.17-------------------------------------------------------------------------#
N <- nrow(data_chap2)
all_store <- factor(c(rep(1, N), rep(2, N), rep(3,N)))
#-2.18-------------------------------------------------------------------------#
TukeyHSD(aov(all_sales ~ all_store))
#-2.19-------------------------------------------------------------------------#
result_1 <- aov(all_sales ~ all_store)
summary(result_1)
#-2.20-------------------------------------------------------------------------#
all_advertising <- rep(data_chap2$Advertising, 3)
#-2.21-------------------------------------------------------------------------#
# 店舗・広告の条件に分けて平均と標準偏差を計算します。
y_mean <- tapply(all_sales,
list(all_advertising,all_store),mean)
y_sd <- tapply(all_sales,
list(all_advertising,all_store),sd)
# まずは棒グラフを描画するオブジェクトbpを定義します。
bp <- barplot(y_mean, beside = TRUE,
ylim = c(0, max(y_mean+y_sd) * 1.1),
col = c("white","grey"),
names = c("Store A", "Store B", "Store C"))
# エラーバー付きの棒グラフを描画します。
arrows(bp, y_mean + y_sd , bp, y_mean - y_sd,
angle = 90, code = 3, length = 0.1)
# 凡例を描画します。
legend("topleft", c("No Ads", "Ads"),fill = c("white", "grey"))
#-2.22-------------------------------------------------------------------------#
result_2 <- aov(all_sales ~ all_store * all_advertising)
summary(result_2)
#------------------------------------------------------------------------------#この章では、パッケージstargazerとlmtestを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
install.packages("stargazer")
install.packages("lmtest")なお、関数install.packagesにオプションとしてrepos = “https://cloud.r-project.org”を付けると、パッケージのファイルはcloudからダウンロードされます。
本文で解説しているコードは以下にあります。開いて確認してください。
#-3.1--------------------------------------------------------------------------#
data_chap3 <- read.csv("chapter_03.csv", header = TRUE)
head(data_chap3, 10)
#-3.2--------------------------------------------------------------------------#
plot(data_chap3$Sales / 1000,
xlab = "time",
ylab = "Sales (1000JPY)",
type = "l")
#-3.3--------------------------------------------------------------------------#
plot(data_chap3$Price, data_chap3$Sales / 1000,
xlab = "Price", ylab = "Sales (1000JPY)", pch = 20)
#-3.4--------------------------------------------------------------------------#
boxplot(Sales / 1000 ~ Display + Feature,
names = c(
"No Feature & No Display", "Display Only",
"Feature Only", "Feature & Display"),
ylab = "Sales (1000JPY)", data = data_chap3)
#-3.5--------------------------------------------------------------------------#
model1 <- lm(Sales ~ Price, data = data_chap3)
model2 <- lm(Sales ~ Price + Display + Feature,
data = data_chap3)
model3 <- lm(Sales ~ Price + Display + Feature + DisFeat,
data = data_chap3)
#-3.6--------------------------------------------------------------------------#
library(stargazer)
stargazer(model1, model2, model3, type = "text")
#-3.7--------------------------------------------------------------------------#
AIC(model1, model2, model3)
BIC(model1, model2, model3)
#-3.8--------------------------------------------------------------------------#
library(lmtest)
lrtest(model1, model2)
lrtest(model2, model3)
#------------------------------------------------------------------------------#この章では、パッケージstargazerを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
install.packages("stargazer")また、以下の関数を定義してください。
#-4.A1-------------------------------------------------------------------------#
roc <- function(y, p){
yo <- y[order(-p)]
n <- length(yo)
curve <- matrix(0, n, 2)
nn <- sum(yo == 0)
pn <- sum(yo == 1)
for (i in 1:n){
curve[i,1] <- sum(yo[1:i] == 0) / nn
curve[i,2] <- sum(yo[1:i] == 1) / pn
}
lst <- as.list(NULL)
lst$curve <- curve
lst$score <- sum(curve[yo == 0,2]) / sum(yo == 0)
return(lst)
}
#-4.A2-------------------------------------------------------------------------#
gain <- function(y, p){
yo <- y[order(-p)]
n <- length(yo)
yg <- yo[1]
for (i in 2:n){
yg[i] <- yg[i - 1] + yo[i]
}
return(yg)
}
#------------------------------------------------------------------------------#本文で解説しているコードは以下にあります。開いて確認してください。
#-4.1--------------------------------------------------------------------------#
#データのインポート
data_chap4 <- read.csv("chapter_04.csv", header = TRUE)
#データ概要の表示
head(data_chap4)
#-4.2--------------------------------------------------------------------------#
# 全体: Aのシェア
mean(data_chap4$YA)
# 全体: Bのシェア
mean(data_chap4$YB)
#-4.3--------------------------------------------------------------------------#
# Aの価格がBの価格より高いとき
mean(data_chap4$YA[data_chap4$PriceA > data_chap4$PriceB])
# Aの価格がBの価格より安いとき
mean(data_chap4$YA[data_chap4$PriceA < data_chap4$PriceB])
#-4.4--------------------------------------------------------------------------#
# Aが特別展示をしているとき
mean(data_chap4$YA[data_chap4$DispA==1])
# Bが特別展示をしているとき
mean(data_chap4$YA[data_chap4$DispB==1])
#-4.5--------------------------------------------------------------------------#
# 訂正前（教科書に掲載）
#data_chap4["PRatio"] <- log(data_chap4[,5]) - log(data_chap4[,4])
# 訂正後
data_chap4["PRatio"] <- log(data_chap4[,4]) - log(data_chap4[,5])
#-4.6--------------------------------------------------------------------------#
formula_1 <- as.formula(YA ~ PRatio + DispA + DispB)
#-4.7--------------------------------------------------------------------------#
model_logit <- glm(formula_1, data = data_chap4,
family = binomial(link = logit))
model_probit <- glm(formula_1, data = data_chap4,
family = binomial(link = probit))
#-4.8--------------------------------------------------------------------------#
summary(model_logit)
summary(model_probit)
#-4.9--------------------------------------------------------------------------#
library(stargazer)
stargazer(model_logit, model_probit, type = "text")
#-4.10-------------------------------------------------------------------------#
formula_p <- as.formula(YA ~ PRatio)
model_logitp <- glm(formula_p, data = data_chap4,
family = binomial(link = logit))
model_probitp <- glm(formula_p, data = data_chap4,
family = binomial(link = probit))
#-4.11-------------------------------------------------------------------------#
AIC(model_logitp, model_logit)
BIC(model_logitp, model_logit)
#-4.12-------------------------------------------------------------------------#
formula_0 <- as.formula(YA ~ 1)
model_logit0 <- glm(formula_0, data = data_chap4,
family = binomial(link = logit))
#-4.13-------------------------------------------------------------------------#
mcfadden_r <- 1 - logLik(model_logit)[1] / logLik(model_logit0)[1]
mcfadden_r
#-4.14-------------------------------------------------------------------------#
data_chap4_c <- data_chap4[data_chap4$Forecast == 0,]
data_chap4_f <- data_chap4[data_chap4$Forecast == 1,]
#-4.15-------------------------------------------------------------------------#
model_logit_c <- glm(formula_1, data = data_chap4_c,
family = binomial(link = logit))
beta <- model_logit_c$coef
#-4.16-------------------------------------------------------------------------#
est_v <- cbind(1,data_chap4_f$PRatio,
data_chap4_f$DispA,
data_chap4_f$DispB) %*% beta
#-4.17-------------------------------------------------------------------------#
est_p <- exp(est_v) / (1 + exp(est_v))
# プロビットモデルなら、est_p <- pnorm(est_v)
#-4.18-------------------------------------------------------------------------#
s <- 0.5 # 閾値が0.5のとき
#s <- mean(data_chap4_c$YA) # 閾値がブランドのシェアのとき
ya_f <- data_chap4_f$YA
ya_ff <- factor(ya_f, levels = c(0, 1))
est_y <- as.numeric(est_p > s)
est_y <- factor(est_y, levels = c(0,1))
m <- table(ya_ff, est_y)
sum(diag(m)) / sum(m)
#-4.19-------------------------------------------------------------------------#
res_roc <- roc(ya_f, est_p)
plot(res_roc$curve, type = "l",
xlab = "False Positive Rate", ylab = "True Positive Rate")
lines(c(0, 1), c(0, 1))
#-4.20-------------------------------------------------------------------------#
res_gain_a <- gain(ya_f, est_p)
plot(c(0, res_gain_a), type = "l",
xlab = "Number of Customers", ylab = "Cumulative Purchase")
lines(c(1, length(ya_f)), c(0, sum(ya_f)))
# 完全予測の点線を追加
res_gain_ap <- gain(ya_f, ya_f)
par(new = TRUE)
plot(c(0, res_gain_ap), type = "l", lty = 2,
ann = FALSE, axes = FALSE)
#-4.21-------------------------------------------------------------------------#
res_gain_b <- gain(1 - ya_f, 1 - est_p)
plot(c(0,res_gain_b), type="l",
xlab = "Number of Customers", ylab = "Cumulative Purchase")
lines(c(1, length(ya_f)), c(0, sum(1 - ya_f)))
# 完全予測の点線を追加
res_gain_ap <- gain(1 - ya_f, 1 - ya_f)
par(new = TRUE)
plot(c(0, res_gain_ap), type = "l", lty = 2,
ann = FALSE, axes = FALSE)この章では、パッケージmlogitを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
install.packages("mlogit")本文で解説しているコードは以下にあります。開いて確認してください。
#-5.1--------------------------------------------------------------------------#
data_chap5 <- read.csv("chapter_05.csv", header = TRUE)
#-5.2--------------------------------------------------------------------------#
barplot(colMeans(data_chap5[,4:7]))
#-5.3--------------------------------------------------------------------------#
barplot(colMeans(data_chap5[,8:11]), main = "PRICE")
#-5.4--------------------------------------------------------------------------#
library("mlogit")
#-5.5--------------------------------------------------------------------------#
data_mlogit <- mlogit.data(data_chap5, varying = 4:23,
shape = "wide")
#-5.6--------------------------------------------------------------------------#
head(data_mlogit)
#-5.7--------------------------------------------------------------------------#
data_est <- data_mlogit[data_mlogit[,3] == 0,]
data_val <- data_mlogit[data_mlogit[,3] == 1,]
#-5.8--------------------------------------------------------------------------#
result_1 <- mlogit(y ~ PRICE + DISPL + FEAT + FD| 1| 0, data_est,
reflevel = "NABISCO")
summary(result_1)
#-5.9--------------------------------------------------------------------------#
result_2 <- mlogit(y ~ PRICE| 1| DISPL + FEAT + FD, data_est,
reflevel = "NABISCO")
summary(result_2)
#-5.10-------------------------------------------------------------------------#
forecast_p <- predict(result_1, data_val)
#-5.11-------------------------------------------------------------------------#
forecast_y <- apply(as.matrix(forecast_p), 1, which.max)
forecast_y <- factor(colnames(forecast_p)[forecast_y],
levels = colnames(forecast_p))
#-5.12-------------------------------------------------------------------------#
observed_y <- factor(data_val[data_val[,5] == 1,4],
levels = colnames(forecast_p))
table(observed_y, forecast_y)
#-5.13-------------------------------------------------------------------------#
sum(diag(table(forecast_y, observed_y))) / length(observed_y)
#------------------------------------------------------------------------------#この章では、パッケージmlogit, flexmixを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
install.packages("mlogit")
install.packages("flexmix")本文で解説しているコードは以下にあります。開いて確認してください。
注）この章で使う関数flexmixによる潜在クラスモデルのパラメータ推定には、PCの性能にも依存しますが、数分～十数分の時間がかかります。その間Rは入力を受け付けないので、十分に時間があるときに実行してください。
#-6.1--------------------------------------------------------------------------#
data_chap6 <- read.csv("chapter_06.csv", header = TRUE)
#-6.2--------------------------------------------------------------------------#
library("mlogit")
cr_data <- mlogit.data(data_chap6,
choice = "choice", id = "id",
varying = 2:13, shape = "wide", sep = ".")
#-6.3--------------------------------------------------------------------------#
head(cr_data)
#-6.4--------------------------------------------------------------------------#
result_logit <- mlogit(choice ~ disp + feat + price,
data = cr_data, reflevel = "private")
summary(result_logit)
#-6.5--------------------------------------------------------------------------#
result_logit_m <- mlogit(choice ~ disp + feat + price,
data = subset(cr_data,sex == 0),
reflevel = "private")
as.data.frame(result_logit_m$coefficients)
result_logit_f <- mlogit(choice ~ disp + feat + price ,
data = subset(cr_data,sex == 1),
reflevel = "private")
as.data.frame(result_logit_f$coefficients)
#-6.6--------------------------------------------------------------------------#
library(flexmix)
result_lc1 <- flexmix(choice ~ disp + feat + price + alt | id,
model = FLXMRcondlogit(strata =~ chid),
data = cr_data, k = 1)
set.seed(123)
result_lc2 <- flexmix(choice ~ disp + feat + price + alt | id,
model = FLXMRcondlogit(strata =~ chid),
data = cr_data, k = 2)
set.seed(123)
result_lc3 <- flexmix(choice ~ disp + feat + price + alt | id,
model = FLXMRcondlogit(strata =~ chid),
data = cr_data, k = 3)
set.seed(123)
result_lc4 <- flexmix(choice ~ disp + feat + price + alt | id,
model = FLXMRcondlogit(strata =~ chid),
data = cr_data, k = 4)
set.seed(123)
result_lc5 <- flexmix(choice ~ disp + feat + price + alt | id,
model = FLXMRcondlogit(strata =~ chid),
data = cr_data, k = 5)
set.seed(123)
result_lc6 <- flexmix(choice ~ disp + feat + price + alt | id,
model = FLXMRcondlogit(strata =~ chid),
data = cr_data, k = 6)
#-6.7--------------------------------------------------------------------------#
BIC(result_lc1, result_lc2, result_lc3,
result_lc4, result_lc5, result_lc6)
#-6.8--------------------------------------------------------------------------#
result_refit_lc4 <- refit(result_lc4)
summary(result_refit_lc4)
#-6.9--------------------------------------------------------------------------#
cid <- cr_data$id
post_lc4 <- unique(cbind(cid, posterior(result_lc4)))
head(post_lc4, 10)
#-6.10-------------------------------------------------------------------------#
result_lc4a <- flexmix(choice ~ disp + feat + price + alt | id,
model = FLXMRcondlogit(strata =~ chid),
concomitant = FLXPmultinom(~ sex + age + fsize),
data = cr_data, k = 4)
parameters(result_lc4a)
#-6.11-------------------------------------------------------------------------#
result_refit_lc4a <- refit(result_lc4a)
summary(result_refit_lc4a, which = "concomitant")
#-6.12-------------------------------------------------------------------------#
plot(result_refit_lc4a)
#------------------------------------------------------------------------------#■ データのダウンロード(chapter_06_choice_data.csv)
chapter_06_choice_data.csvの概要
| 項目 | 概要 |
|---|---|
| Cid | 顧客のID番号 |
| feat.yoplait | ブランドyoplaitのチラシ広告の有無を示すダミー変数 |
| feat.dannon | ブランドyoplaitのチラシ広告の有無を示すダミー変数 |
| feat.hiland | ブランドyoplaitのチラシ広告の有無を示すダミー変数 |
| feat.weight | ブランドyoplaitのチラシ広告の有無を示すダミー変数 |
| price.yoplait | ブランドyoplaitの価格 |
| price.dannon | ブランドyoplaitの価格 |
| price.hiland | ブランドyoplaitの価格 |
| price.weight | ブランドyoplaitの価格 |
| choice | 選択されたブランド |
| sex | 顧客の性別を示すダミー変数（０＝女性，１＝男性） |
| age | 顧客の年齢 |
| married | 顧客の婚姻状況（０＝結婚していない，１＝結婚している） |
この章では、パッケージcensReg, sampleSelectionを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
install.packages(“censReg”) install.packages(“sampleSelection”)
本文で解説しているコードは以下にあります。開いて確認してください。
<details>
<summary>第7章のコードはこちら</summary>
```r
#-7.1--------------------------------------------------------------------------#
data_chap7 <- read.csv("chapter_07.csv", header = TRUE)
head(data_chap7, 10)
#-7.2--------------------------------------------------------------------------#
table(data_chap7$Purchase)
#-7.3--------------------------------------------------------------------------#
hist(data_chap7$Amount, xlab="Amount",
ylab = "Number of Customers", main = "")
#-7.4--------------------------------------------------------------------------#
mean(data_chap7$Amount[data_chap7$Amount > 0])
median(data_chap7$Amount[data_chap7$Amount > 0])
#-7.5--------------------------------------------------------------------------#
opar <- par(no.readonly = TRUE)
layout(matrix(c(1:6), 2, 3 , byrow = "TRUE"))
plot(data_chap7$Sex, data_chap7$Amount,
xlab = "Sex", ylab = "Amount")
plot(data_chap7$Age, data_chap7$Amount,
xlab = "Age", ylab = "Amount")
plot(data_chap7$Fsize, data_chap7$Amount,
xlab = "Fsize", ylab = "Amount")
plot(data_chap7$Income, data_chap7$Amount,
xlab = "Income", ylab = "Amount")
plot(data_chap7$Ownhouse, data_chap7$Amount,
xlab = "Ownhouse", ylab = "Amount")
par(opar)
#-7.6--------------------------------------------------------------------------#
## 全データを用いた線形回帰モデル
result_reg1 <- lm(Amount ~ Sex + Age + Fsize + Income +
Ownhouse + Crossbuying,
data = data_chap7)
summary(result_reg1)
## オンライン店舗の利用者のデータを用いた線形回帰モデル
result_reg2 <- lm(Amount ~ Sex + Age + Fsize + Income +
Ownhouse + Crossbuying,
data = subset(data_chap7, Amount > 0))
summary(result_reg2)
## タイプIトービットモデル
library(censReg)
result_tobit1 <- censReg(Amount ~ Sex + Age + Fsize + Income +
Ownhouse + Crossbuying,
data = data_chap7)
summary(result_tobit1)
#-7.7--------------------------------------------------------------------------#
margEff(result_tobit1)
#-7.8--------------------------------------------------------------------------#
library(sampleSelection)
result_tobit2 <- heckit(
selection = Purchase ~ Sex + Age + Crossbuying + Pfreq + Rduration,
outcome = Amount ~ Sex + Age + Fsize + Income + Ownhouse,
data = data_chap7, method = "ml")
summary(result_tobit2)
#-7.9--------------------------------------------------------------------------#
#トービットモデルの推定値
res_coeff <- coef(result_tobit2)
#選択モデルの推定値
res_theta <- as.vector(res_coeff[1:6])
z_data <- as.matrix(
cbind(Intercept = 1,
subset(data_chap7,
select = c(Sex, Age, Crossbuying, Pfreq, Rduration)
)
)
)
z_theta <- crossprod(t(z_data), res_theta)
#逆ミルズ比
inv_mills_r <- dnorm(z_theta) / pnorm(z_theta)
#選択モデルでの年齢の係数
theta_age <- res_coeff[3]
#購買金額モデルでの年齢の係数
beta_age<-res_coeff[9]
res_sigma<-res_coeff[13]
res_rho<-res_coeff[14]
#年齢の限界効果
marg_effect <- beta_age - res_rho * res_sigma * theta_age *
inv_mills_r * (z_theta+inv_mills_r)
mean(marg_effect)
#------------------------------------------------------------------------------#■ データのダウンロード(data/chapter_07_loyalty.csv)
chapter_07_loyalty.csvの概要
| 項目 | 概要 |
|---|---|
| Cid | 顧客のID番号 |
| Loyalty | ロイヤルティ・プログラムへの参加を示すダミー変数（０＝参加していない，１＝参加している） |
| Pamount | 購買金額の対数 |
| Sex | 顧客の性別を示すダミー変数（０＝女性，１＝男性） |
| Age | 顧客の年齢 |
| Fsize | 顧客の家族人数 |
| Logincome | 顧客の収入の対数 |
この章では、パッケージpsclを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
install.packages("pscl")本文で解説しているコードは以下にあります。開いて確認してください。
#-8.1--------------------------------------------------------------------------#
data_chap8 <- read.csv("chapter_08.csv", header = TRUE)
head(data_chap8)
#-8.2--------------------------------------------------------------------------#
hist(data_chap8$Contracts, 20)
#-8.3--------------------------------------------------------------------------#
y <- data_chap8$Contracts
x2 <- data_chap8$Discount > 0
t.test(y[x2], y[!x2])
#-8.4--------------------------------------------------------------------------#
y <- data_chap8$Contracts
x1 <- data_chap8$Advertising > 0
t.test(y[x1], y[!x1])
#-8.5--------------------------------------------------------------------------#
result_poi0 <- glm(Contracts ~ Advertising + Discount +
Holiday + RainTotal + TempAve,
family = poisson(link = log),
data = data_chap8)
#-8.6--------------------------------------------------------------------------#
summary(result_poi0)
#-8.7--------------------------------------------------------------------------#
HLTemp <- abs(data_chap8$TempAve - mean(data_chap8$TempAve))
data_chap8["HLTemp"] <- HLTemp
#-8.8--------------------------------------------------------------------------#
result_poi1 <- glm(Contracts ~ Advertising + Discount +
Holiday + RainTotal + HLTemp,
family = poisson(link = log),
data = data_chap8)
#-8.9--------------------------------------------------------------------------#
AIC(result_poi0, result_poi1)
#-8.10-------------------------------------------------------------------------#
summary(result_poi1)
#-8.11-------------------------------------------------------------------------#
library(MASS)
result_nb <- glm.nb(Contracts ~ Advertising + Discount +
Holiday + RainTotal + HLTemp,
link=log,
data = data_chap8)
#-8.12-------------------------------------------------------------------------#
summary(result_nb)
#-8.13-------------------------------------------------------------------------#
AIC(result_poi1, result_nb)
#-8.14-------------------------------------------------------------------------#
table(data_chap8$Contracts > 0)
mean(data_chap8$Contracts > 0)
#-8.15-------------------------------------------------------------------------#
library(pscl)
#-8.16-------------------------------------------------------------------------#
result_zip <- zeroinfl(Contracts ~
Advertising + Discount + Holiday + RainTotal + HLTemp |
Advertising + Discount + Holiday + RainTotal + HLTemp,
dist = "poisson", data = data_chap8)
#-8.17-------------------------------------------------------------------------#
result_zinb <- zeroinfl(Contracts ~
Advertising + Discount + Holiday + RainTotal + HLTemp |
Advertising + Discount + Holiday + RainTotal + HLTemp,
dist = "negbin", data = data_chap8)
#-8.18-------------------------------------------------------------------------#
AIC(result_zip, result_zinb)
#-8.19-------------------------------------------------------------------------#
summary(result_zinb)
#------------------------------------------------------------------------------#この章では、パッケージsurvival, stargazarを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
install.packages("survival")
install.packages("stargazar")本文で解説しているコードは以下にあります。開いて確認してください。
#-9.1--------------------------------------------------------------------------#
data_chap9 <- read.csv("chapter_09.csv", header = TRUE)
head(data_chap9)
#-9.2--------------------------------------------------------------------------#
nrow(subset(data_chap9, Event == 0))
#-9.3--------------------------------------------------------------------------#
summary(data_chap9[data_chap9$Event == 1,]$Duration)
#-9.4--------------------------------------------------------------------------#
hist(data_chap9$Duration, main = "",
xlab = "Duration", ylab = "Number of Cars")
#-9.5--------------------------------------------------------------------------#
dur_ads <- data_chap9[data_chap9$Ads == 1,1]
dur_no_ads <- data_chap9[data_chap9$Ads == 0,1]
boxplot(dur_ads, dur_no_ads,
names = c("Ads", "No Ads"), ylab = "Duration")
#-9.6--------------------------------------------------------------------------#
library(survival)
result_haz <- survfit(Surv(Duration, Event) ~ Ads, data = data_chap9)
result_haz
#-9.7--------------------------------------------------------------------------#
summary(result_haz)
#-9.8--------------------------------------------------------------------------#
plot(result_haz, lty = c(1, 2),
xlab = "Duration", ylab = "Survival Function")
legend("topright", c("No Ads", "Ads"), lty = c(1, 2))
#-9.9--------------------------------------------------------------------------#
survdiff(Surv(Duration, Event) ~ Ads, data = data_chap9)
#-9.10-------------------------------------------------------------------------#
result_cox <- coxph(Surv(Duration,Event) ~ Ads + Price + Carage +
Origin + Trans + Mileage, data = data_chap9)
summary(result_cox)
#-9.11-------------------------------------------------------------------------#
# Exponential model
result_ex <- survreg(Surv(Duration, Event) ~ Ads + Price + Carage +
Origin + Trans + Mileage, data = data_chap9,
dist = "exponential")
# Weibull model
result_wei <- survreg(Surv(Duration, Event) ~ Ads + Price + Carage +
Origin + Trans + Mileage, data = data_chap9,
dist = "weibull")
# Logistic model
result_log <- survreg(Surv(Duration, Event) ~ Ads + Price + Carage +
Origin + Trans + Mileage, data = data_chap9,
dist = "logistic")
# Lognormal model
result_ln <- survreg(Surv(Duration, Event) ~ Ads + Price + Carage +
Origin + Trans + Mileage, data = data_chap9,
dist = "lognormal")
# Loglogistic model
result_ll <- survreg(Surv(Duration, Event) ~ Ads + Price + Carage +
Origin + Trans + Mileage, data = data_chap9,
dist = "loglogistic")
#-9.12-------------------------------------------------------------------------#
library(stargazer)
stargazer(result_ex, result_wei, result_log, result_ln, result_ll,
type = "text")
#-9.13-------------------------------------------------------------------------#
BIC(result_ex, result_wei, result_log, result_ln, result_ll)
#-9.14-------------------------------------------------------------------------#
predict_ln <- predict(result_ln, newdata = data_chap9)
plot(data_chap9$Duration, predict_ln,
xlab = "Observation", ylab = "Prediction")
#------------------------------------------------------------------------------#■ データのダウンロード(chapter_09_churn_data.csv)
chapter_09_churn_data.csvの概要
| 項目 | 概要 |
|---|---|
| Cid | 顧客のID番号 |
| Duration | 解約までの期間（週数） |
| Event | 解約の有無を示すダミー変数（０＝解約なし，１＝解約あり） |
| DM | 顧客がダイレクトメール（DM）を受け取っているかどうかのダミー変数（1=受け取っている，0=受け取っていない） |
| Sex | 顧客の性別を示すダミー変数（０＝女性，１＝男性） |
| Age | 顧客の年齢 |
| Freq | 顧客のサービス利用頻度 |
| Amount | 一か月あたりのサービス利用金額 |
■ データのダウンロード(chapter_10_1.csv)
■ データのダウンロード(chapter_10_2.csv)
この章では、パッケージconjointを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
install.packages("conjoint")本文で解説しているコードは以下にあります。開いて確認してください。
#-10.1-------------------------------------------------------------------------#
data_chap10_1 <- read.csv("chapter_10_1.csv", header = TRUE)
head(data_chap10_1)
#-10.2-------------------------------------------------------------------------#
result_lm <- lm( PI ~ Age + Sex + Sat + Qual, data = data_chap10_1)
summary(result_lm)
#-10.3-------------------------------------------------------------------------#
plot(data_chap10_1$Age, data_chap10_1$PI, xlab = "Age",
ylab = "Repeat Purchase", pch = 16)
abline(lm(PI ~ Age, data = data_chap10_1))
#-10.4-------------------------------------------------------------------------#
data_chap10_1$PI <- as.ordered(data_chap10_1$PI)
head(data_chap10_1$PI)
#-10.5-------------------------------------------------------------------------#
library(MASS)
result_ol <- polr( PI ~ Age + Sex + Sat + Qual, data = data_chap10_1)
summary(result_ol)
#-10.6-------------------------------------------------------------------------#
proby <- fitted.values(result_ol)
head(proby)
#-10.7-------------------------------------------------------------------------#
library(conjoint)
set.seed(123)
caratt <- list(price = c("2M JPY", "2.5M JPY"),
fuel = c("10㎞/l", "14㎞/l"),
displacement = c("1500㏄", "2000㏄"))
caratt
#-10.8-------------------------------------------------------------------------#
all_com <- expand.grid(caratt)
all_com
#-10.9-------------------------------------------------------------------------#
ex_design <- caFactorialDesign(data = all_com, type = "orthogonal")
ex_design
#-10.10------------------------------------------------------------------------#
data_chap10_2 <- read.csv("chapter_10_2.csv", header = TRUE)
head(data_chap10_2)
#-10.11------------------------------------------------------------------------#
car_data <- NULL
for(i in 1:nrow(data_chap10_2)){
pref_order <- as.vector(t(data_chap10_2[i,2:5]))
subdat <- cbind(pref_order, ex_design)
car_data <- rbind(car_data, subdat)
}
head(car_data)
#-10.12------------------------------------------------------------------------#
car_data$pref_order <- as.ordered(car_data$pref_order)
head(car_data$pref_order)
#-10.13------------------------------------------------------------------------#
result_car <- polr(pref_order ~ price + fuel + displacement,
data = car_data)
summary(result_car, digits = 3)
#-10.14------------------------------------------------------------------------#
car_range <- abs(result_car$coefficients[1:3])
car_range
#-10.15------------------------------------------------------------------------#
importance <- car_range / sum(car_range)
importance
#------------------------------------------------------------------------------#■ データのダウンロード(chapter_10_conjoint_data.csv)
chapter_10_conjoint_data.csvの概要
| 項目 | 概要 |
|---|---|
| 順位 | 選好の順位 |
| 価格 | PCの価格 |
| メモリ | PCのメモリ |
| 重量 | PCの重量 |
この章では、パッケージlavaanを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
install.packages("lavaan")本文で解説しているコードは以下にあります。開いて確認してください。
#-11.1-------------------------------------------------------------------------#
data_chap11 <- read.csv("chapter_11.csv", header = TRUE)
#-11.2-------------------------------------------------------------------------#
dat_lov <- data_chap11[,1:9]
#-11.3-------------------------------------------------------------------------#
cor(dat_lov)
#-11.4-------------------------------------------------------------------------#
eigenvalues <- eigen(cor(dat_lov))$values
plot(eigenvalues, type = "b")
abline(h = 1)
#-11.5-------------------------------------------------------------------------#
result_fa <- factanal(dat_lov, 2)
result_fa <- factanal(dat_lov, 2)
#-11.6-------------------------------------------------------------------------#
label_lov <- c("achievement factor", "belongingness factor")
fa_loadings <- result_fa$loadings[,]
colnames(fa_loadings) <- label_lov
plot(fa_loadings, pch = "")
text(fa_loadings[,1], fa_loadings[,2], row.names(fa_loadings))
#-11.7-------------------------------------------------------------------------#
result_fa_n <- factanal(dat_lov, 2, rotation = "none")
fa_loadings_n <- result_fa_n$loadings[,]
plot(fa_loadings_n, pch = "")
text(fa_loadings_n[,1], fa_loadings_n[,2], row.names(fa_loadings_n))
#-11.8-------------------------------------------------------------------------#
result_fa <- factanal(dat_lov, 2, scores = "regression")
fa_scores <- result_fa$scores
colnames(fa_scores) <- label_lov
plot(fa_scores)
#-11.9-------------------------------------------------------------------------#
dat_lu <- data_chap11[,10:17]
#-11.11------------------------------------------------------------------------#
cor(dat_lu)
plot(eigen(cor(dat_lu))$values, type = "b")
abline(h = 1)
#-11.12------------------------------------------------------------------------#
result_fa <- factanal(dat_lu, 2)
fa_loadings <- result_fa$loadings[,]
plot(fa_loadings, pch = "")
text(fa_loadings[,1], fa_loadings[,2], row.names(fa_loadings))
#-11.13------------------------------------------------------------------------#
library("lavaan")
#-11.14------------------------------------------------------------------------#
model_cfa <- '
AT =~ LUAT1 + LUAT2 + LUAT3
BE =~ LUBE1 + LUBE2 + LUBE3 + LUBE4 + LUBE5
'
#-11.15------------------------------------------------------------------------#
result_cfa <- lavaan::sem(model_cfa, data = dat_lu)
fitMeasures(result_cfa)
summary(result_cfa, standardized = TRUE)
#-11.16------------------------------------------------------------------------#
fitMeasures(result_cfa)
#-11.17------------------------------------------------------------------------#
summary(result_cfa, standardized = TRUE)
#-11.18------------------------------------------------------------------------#
alpha.coef <- function(x){
z <- (ncol(x) / (ncol(x) - 1)) *
(1 - sum(diag(var(x))) / var(rowSums(x)))
return(z)
}
#-11.19------------------------------------------------------------------------#
alpha.coef(dat_lu[,1:3])
alpha.coef(dat_lu[,4:8])
#-11.20------------------------------------------------------------------------#
standardizedSolution(result_cfa)
#-11.21------------------------------------------------------------------------#
cr_ave <- function(Flab, StdEst){
FL <- StdEst[(StdEst["lhs"] == Flab) & (StdEst["op"] == "=~"),3]
FLlist <- colSums(apply(StdEst["rhs"], 1,
function(x,fl){x==fl},FL)) > 0
Lm <- StdEst[FLlist & (StdEst["op"] == "=~"), 4]
e <- StdEst[FLlist & (StdEst["op"] == "~~"), 4]
CR <- sum(Lm)^2 / (sum(Lm)^2 + sum(e))
AVE <- sum(Lm^2) / (sum(Lm^2) + sum(e))
return(c(CR,AVE))
}
#-11.22------------------------------------------------------------------------#
std_est <- standardizedSolution(result_cfa)
cr_ave("AT", std_est)
cr_ave("BE", std_est)
#-11.23------------------------------------------------------------------------#
standardizedSolution(result_cfa)
#------------------------------------------------------------------------------#注）第12章では、第11章と同じデータを使います。
この章では、パッケージlavaanを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
install.packages("lavaan")本文で解説しているコードは以下にあります。開いて確認してください。
#-12.1-------------------------------------------------------------------------#
library("lavaan")
#-12.2-------------------------------------------------------------------------#
data_chap12 <- read.csv(file="chapter_11.csv", header = TRUE)
data_lm <- data.frame(Z = data_chap12[,18],
Y = data_chap12[,19], X = data_chap12[,20])
#-12.3-------------------------------------------------------------------------#
model_lm <- 'Y ~ 1 + X + Z'
#-12.4-------------------------------------------------------------------------#
result_lm <- lavaan::sem(model_lm, data = data_lm)
#-12.5-------------------------------------------------------------------------#
summary(result_lm, standardized = FALSE)
#-12.6-------------------------------------------------------------------------#
summary(lm(Y ~ X + Z, data = data_lm))
#-12.7-------------------------------------------------------------------------#
model_sur <- 'X ~ 1 + Z
Y ~ 1 + Z
X ~~ Y'
#-12.8-------------------------------------------------------------------------#
result_sur <- lavaan::sem(model_sur, data = data_lm)
summary(result_sur, standardized = FALSE)
#-12.9-------------------------------------------------------------------------#
model_ie <- 'X ~ 1 + a12 * Y
Y ~ 1 + a23 * Z
a12a13 := a12 * a23'
#-12.10------------------------------------------------------------------------#
result_ie <- lavaan::sem(model_ie, data = data_lm)
summary(result_ie, standardized = TRUE)
#-12.11------------------------------------------------------------------------#
result_iebs <- lavaan::sem(model_ie, data = data_lm,
se = "bootstrap", bootstrap = 2000)
summary(result_iebs, standardized = TRUE)
#-12.11------------------------------------------------------------------------#
model_ie2 <- 'X ~ 1 + a12 * Y + a13 * Z
Y ~ 1 + a23 * Z
aIndirect := a12 * a23
aTotal := a12 * a23 + a13'
#-12.12------------------------------------------------------------------------#
result_ie2 <- lavaan::sem(model_ie2, data = data_lm)
summary(result_ie2, standardized = TRUE)
#-12.13------------------------------------------------------------------------#
model_2f <- 'AT =~ LUAT1 + LUAT2 + LUAT3
BE =~ LUBE1 + LUBE2 + LUBE3 + LUBE5
LU =~ AT + BE
LU ~~ 1 * LU'
#-12.14------------------------------------------------------------------------#
result_2f <- lavaan::sem(model_2f, data = data_chap12)
summary(result_2f, standardized = TRUE)
#-12.15------------------------------------------------------------------------#
standardizedSolution(result_2f)
#-12.16------------------------------------------------------------------------#
model_nfc <- 'AT =~ LUAT1 + LUAT2 + LUAT3
BE =~ LUBE1 + LUBE2 + LUBE3 + LUBE5
NFC =~ NFC1 + NFC2 + NFC3 + NFC4
AT ~ NFC
BE ~ NFC'
#-12.17------------------------------------------------------------------------#
result_nfc <- lavaan::sem(model_nfc, data = data_chap12)
summary(result_nfc, standardized = TRUE)
#-12.18------------------------------------------------------------------------#
result_nfc_mg <- lavaan::sem(model_nfc, group = "gen",
data = data_chap12)
#-12.19------------------------------------------------------------------------#
model_nfc2 <- 'AT =~ LUAT1 + LUAT2 + LUAT3
BE =~ LUBE1 + LUBE2 + LUBE3 + LUBE5
NFC =~ NFC1 + NFC2 + NFC3 + NFC4
AT ~ c(a11, a12) * NFC
BE ~ c(a21, a22) * NFC
AT ~~ c(v1, v2) * BE'
#-12.20------------------------------------------------------------------------#
result_nfc_mg2 <- lavaan::sem(model_nfc2, group = "gen", data = data_chap12)
summary(result_nfc_mg2, standardized = TRUE)
#------------------------------------------------------------------------------#この章では、パッケージbaysm, invgamma, rstanarm, bayestestRを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
(注)これらのパッケージは初回のダウンロードとインストールに長い時間を要する可能性がありますので、Google Colabで本章の分析を行うのは推奨しません。
install.packages("bayesm")
install.packages("invgamma")
install.packages("rstanarm")
install.packages("bayestestR")注）パッケージのダウンロードには長い時間がかかる場合があります。
本文で解説しているコードは以下にあります。開いて確認してください。
#-13.1-------------------------------------------------------------------------#
library(bayesm)
data(cheese)
data_cheese <- cheese
data_cheese$VOLUME <- log(data_cheese$VOLUME)
data_cheese$PRICE <- log(data_cheese$PRICE)
head(data_cheese)
#-13.2-------------------------------------------------------------------------#
library(invgamma)
# 事前分布の設定
m_0 <- 8
V_0 <- 20 # 1/sigma^2
a_0 <- 8
b_0 <- 8
# データの設定
set.seed(1)
Y <- rnorm(1000, mean = 10, sd = 2)
y_bar <- mean(Y)
n <- length(Y)
# 事後分布のパラメータ
V_n <- 1 / (1/V_0 + n)
m_n <- V_n * (1/V_0 * m_0 + n*y_bar)
a_n <- a_0 + n/2
b_n <- b_0 + 1/2 * (m_0^2 * 1/V_0 + sum(Y^2) - m_n^2 * 1/V_n)
## MCMCシミュレーションの設定
r <- 2000 # イタレーションの回数
burn <- 1000 # バーンイン期間
ret<- r - burn # 保存する標本数
Mdraw <- rep(0, ret)
Sigmadraw <- rep(0, ret)
for (j in 1:r){
Sigmanew <- rinvgamma(1, shape = a_n, rate = b_n)
Mnew <- rnorm(1, mean = m_n, sd = sqrt(V_n * Sigmanew))
if(j > burn){
Sigmadraw[j - burn] <- Sigmanew
Mdraw[j - burn] <- Mnew
}
}
#-13.3-------------------------------------------------------------------------#
mat <- quantile(Mdraw,probs = c(0.5, 0.025, 0.975))
mat
mat <- quantile(Sigmadraw,probs = c(0.5, 0.025, 0.975))
mat
#-13.4-------------------------------------------------------------------------#
library(rstanarm)
Regmodel <- stan_glm(VOLUME ~ DISP + PRICE, data = data_cheese,
chains = 1, seed = 1234)
#-13.5-------------------------------------------------------------------------#
library(bayestestR)
describe_posterior(Regmodel)
#-13.6--------------------------------------------------------------------------#
data_chap13 <- read.csv("chapter_13.csv", header = TRUE)
head(data_chap13)
#-13.7-------------------------------------------------------------------------#
nreg <- nrow(data_chap13)
Cons <- c(rep(1, nreg))
Hist <- data_chap13$Hist
Loy <- data_chap13$Loy
Z <- cbind(Cons, Hist, Loy)
head(Z)
#-13.8-------------------------------------------------------------------------#
Ret <- unique(data_chap13$Ret)
regdata <- NULL
for (i in 1:nreg) {
subdat <- subset(data_cheese, RETAILER == Ret[i])
Cons <- rep(1,nrow(subdat))
X <- cbind(Cons, subset(subdat, select = c(DISP, PRICE)))
X <- as.matrix(X)
y <- subset(subdat, select = c(VOLUME))
y <- as.matrix(y)
regdata[[i]] <- list(y = y, X = X)
}
Data <- list(regdata = regdata, Z = Z)
#-13.9-------------------------------------------------------------------------#
set.seed(2)
Mcmc <- list(R = 2000, keep = 1)
out <- rhierLinearModel(Data = Data, Mcmc = Mcmc)
#-13.10------------------------------------------------------------------------#
mat <- apply(out$Deltadraw, 2, quantile, probs = c(0.5, 0.025, 0.975))
mat <- t(mat)
colnames(mat) <- c("median", "Lower 2.5%", "Upper 97.5%")
print(mat)
#-13.11------------------------------------------------------------------------#
Delta_estimates <- matrix(mat[,1], 3, 3)
colnames(Delta_estimates) <- c("CONS", "DISP", "PRICE")
rownames(Delta_estimates) <- c("CONS", "HIST", "LOY")
print(Delta_estimates)
#-13.12------------------------------------------------------------------------#
plot(out$betadraw)
#------------------------------------------------------------------------------#■ データのダウンロード(chapter_14_1.csv)
■ データのダウンロード(chapter_14_2.csv)
この章では、パッケージbaysm, rstanarm, bayestestRを使います。最初に以下のコマンドを実行してパッケージをダウンロードしておいてください。
(注)これらのパッケージは初回のダウンロードとインストールに長い時間を要する可能性がありますので、Google Colabで本章の分析を行うのは推奨しません。
install.packages("bayesm")
install.packages("rstanarm")
install.packages("bayestestR")本文で解説しているコードは以下にあります。開いて確認してください。
#-14.1-------------------------------------------------------------------------#
data_chap14_1 <- read.csv("chapter_14_1.csv", header = TRUE)
head(data_chap14_1, 20)
#-14.2-------------------------------------------------------------------------#
plot(data_chap14_1$Pamount, data_chap14_1$Pdur,
xlab = "Pamount", ylab = "Pdur")
#-14.3-------------------------------------------------------------------------#
cor(data_chap14_1$Pamount, data_chap14_1$Pdur)
#-14.4-------------------------------------------------------------------------#
cor.test(data_chap14_1$Pamount,data_chap14_1$Pdur)
#-14.5-------------------------------------------------------------------------#
library(rstanarm)
result_Pamount <- stan_glm(
Pamount ~ Disc + Event + Temp + Sqtemp + Rain,
data = data_chap14_1, chains = 2, seed = 1234)
#-14.6-------------------------------------------------------------------------#
library(bayestestR)
describe_posterior(result_Pamount)
#-14.7-------------------------------------------------------------------------#
result_Pdur <- stan_glm(
Pdur ~ Disc + Event + Temp + Sqtemp + Rain,
data = data_chap14_1, chains = 2, seed = 1234)
describe_posterior(result_Pdur)
#-14.8-------------------------------------------------------------------------#
## 従属変数
Y1 <- data_chap14_1$Pamount
Y2 <- data_chap14_1$Pdur
## 独立変数
Cons <- rep(1, nrow(data_chap14_1))
X <- as.matrix(data.frame(Cons,
data_chap14_1[,c(4:ncol(data_chap14_1))])
)
m <- 2 #従属変数の数
k <- ncol(X) #切片を含めた独立変数の数
#-14.9-------------------------------------------------------------------------#
## 事前分布の設定
betabar <- rep(0, k*m)
Bbar <- matrix(betabar, ncol = m)
A <- diag(rep(0.01, k))
nu <- m
V <- nu * diag(m)
#-14.10------------------------------------------------------------------------#
## MCMCsシミュレーションの設定
r <- 2000 # イタレーションの回数
burn <- 1000 # 除去する最初のサンプル数
ret <- r - burn # 保存するサンプル数
#-14.11------------------------------------------------------------------------#
betadraw <- matrix(double(ret*k*m), ncol = k*m)
Sigmadraw <- matrix(double(ret*m*m), ncol = m*m)
#-14.12------------------------------------------------------------------------#
library(bayesm)
set.seed(1)
for (j in 1:r){
Y <- cbind(Y1, Y2)
# パラメータのサンプリング
out <- rmultireg(Y, X, Bbar, A, nu, V)
## 新しいパラメータ値の保存
if(j > burn){
betadraw[j - burn,] <- out$B
Sigmadraw[j - burn,] <- out$Sigma
}
}
#-14.13------------------------------------------------------------------------#
mat <- apply(betadraw, 2, quantile, probs = c(0.5, 0.025, 0.975))
mat <- t(mat)
#購買金額の推定結果
beta1 <- mat[1:k,]
rownames(beta1) <- colnames(X)
colnames(beta1) <- c("Median", "CI_2.5%", "CI_97.5%")
print(beta1, digits = 3)
#-14.14------------------------------------------------------------------------#
#滞在時間の推定結果
beta2 <- mat[(k+1):(2*k),]
rownames(beta2) <- colnames(X)
colnames(beta2) <- c("Median", "CI_2.5%", "CI_97.5%")
print(beta2, digits = 3)
#-14.15------------------------------------------------------------------------#
#分散共分散の推定結果
mas <- apply(Sigmadraw, 2 ,quantile, probs = c(0.025, 0.5, 0.975))
mas <- t(mas)
print(mas)
covmat <- matrix(mas[,2], m, m)
#相関行列
cormat <- cov2cor(covmat)
print(cormat)
#-14.16------------------------------------------------------------------------#
data_chap14_2 <- read.csv("chapter_14_2.csv", header = TRUE)
head(data_chap14_2)
#-14.17------------------------------------------------------------------------#
# 従属変数の散布図
pairs(data_chap14_2[2:4])
#-14.18------------------------------------------------------------------------#
# 相関行列
cor(data_chap14_2[2:4])
#-14.19------------------------------------------------------------------------#
## 従属変数
Y1 <- data_chap14_2$TV
Y2 <- data_chap14_2$Internet
Y3 <- data_chap14_2$Ecsite
#-14.20------------------------------------------------------------------------#
## 独立変数
Cons <- rep(1, nrow(data_chap14_2))
X1 <- as.matrix(data.frame(Cons, data_chap14_2[,c(5:9)]))
X2 <- as.matrix(data.frame(Cons, data_chap14_2[,c(5:10)]))
X3 <- as.matrix(data.frame(Cons, data_chap14_2[,c(5:11)]))
m <- 3 #従属変数の数
k <- ncol(X1) + ncol(X2) + ncol(X3) #切片を含めた独立変数の数
## 事前分布のパラメータ
betabar <- rep(0, k)
A <- diag(rep(0.01, k))
nu <- m
V <- nu * diag(m)
## MCMCシミュレーションの設定
r <- 2000 # イタレーションの回数
burn <- 1000 # 除去する最初のサンプル数
ret<- r - burn # 保存するサンプル数
betadraw <- matrix(double(ret*k), ncol = k)
Sigmadraw <- matrix(double(ret*m*m), ncol = m*m)
library(bayesm)
set.seed(2)
for (j in 1:r){
regdata <- NULL
regdata[[1]] <- list(y = Y1, X = X1)
regdata[[2]] <- list(y = Y2, X = X2)
regdata[[3]] <- list(y = Y3, X = X3)
# Bayesian SURモデルのサンプリング
out <- rsurGibbs(Data = list(regdata = regdata),
Mcmc = list(R = 1, nprint = 0))
## 新しいパラメータ値の保存
if(j > burn){
betadraw[j - burn,] <- out$betadraw
Sigmadraw[j - burn,] <- out$Sigmadraw
}
}
#-14.21------------------------------------------------------------------------#
mat <- apply(betadraw, 2, quantile, probs = c(0.5, 0.025, 0.975))
mat <- t(mat)
#-14.22------------------------------------------------------------------------#
#テレビ視聴時間の推定結果
beta1 <- mat[1:6,]
rownames(beta1) <- colnames(X1)
colnames(beta1) <- c("Median", "CI_2.5%", "CI_97.5%")
print(beta1, digits = 3)
#-14.23------------------------------------------------------------------------#
#インターネット閲覧時間の推定結果
beta2 <- mat[7:13,]
rownames(beta2) <- colnames(X2)
colnames(beta1) <- c("Median", "CI_2.5%", "CI_97.5%")
print(beta2, digits = 3)
#-14.24------------------------------------------------------------------------#
#ECサイト閲覧時間の推定結果
beta3 <- mat[14:21,]
rownames(beta3) <- colnames(X3)
colnames(beta3) <- c("Median", "CI_2.5%", "CI_97.5%")
print(beta3, digits = 3)
#-14.25------------------------------------------------------------------------#
#分散共分散の推定結果
mas <- apply(Sigmadraw, 2 ,quantile, probs = c(0.025, 0.5, 0.975))
mas <- t(mas)
print(mas)
covmat <- matrix(mas[,2], m, m)
#相関行列
cormat <- cov2cor(covmat)
print(cormat, digits = 3)
#------------------------------------------------------------------------------#■ データのダウンロード(chapter_14_mtvreg.csv)
chapter_14_mtvreg.csvの概要
| 項目 | 概要 |
|---|---|
| Cid | 顧客のID番号 |
| Visit | 顧客の一定期間における来店頻度 |
| Pamount | 来店当たりの購買金額 |
| Sex | 顧客の性別を示すダミー変数（０＝女性，１＝男性） |
| Age | 顧客の年齢 |
| Crossbuying | 来店当たりの購買品数 |
| Rdur | 入会期間 |
■ コーパスデータのダウンロード(chapter_15_corpus.xlsx)
この章では、パッケージopenxlsx,RMeCabを使います。ただし、RMeCabのインストール方法は環境によって異なるので、それぞれの環境に合った方法を使ってください。
まず、openxlsxのインストールはどの環境でも共通です。
install.packages("openxlsx")RMeCabのインストールについては、PCにRをダウンロード・インストールして使う場合とGoogle Colabの場合で異なります。また、PCにインストールしたRを使う場合には、事前にMeCabをインストールしておく必要があります。
#------------------------------------------------------------------------------#
# Windows, Mac OSなどにインストールしたRを使う場合
install.packages("RMeCab", repos = "https://rmecab.jp/R")
#------------------------------------------------------------------------------#
# Google Colabの場合
# MeCabをインストールしてからRMeCabをインストールします
# 参考: https://qiita.com/satoru-ka/items/f47ecd43d9e69d64b0a8
system('sudo apt-get install -y mecab', intern=TRUE)
system('sudo apt-get install -y libmecab-dev', intern=TRUE)
system('sudo apt-get install -y mecab-ipadic-utf8', intern=TRUE)
install.packages("RMeCab", repos = "https://rmecab.jp/R")
#------------------------------------------------------------------------------#本文で解説しているコードは以下にあります。開いて確認してください。
#-15.1-------------------------------------------------------------------------#
library(openxlsx)
data_chap15 <- read.xlsx("chapter_15.xlsx", sheet = "Sheet1")
#-15.2-------------------------------------------------------------------------#
colnames(data_chap15)
#-15.3-------------------------------------------------------------------------#
data_chap15$Content[1]
#-15.4-------------------------------------------------------------------------#
library(RMeCab)
#-15.5-------------------------------------------------------------------------#
RMeCabC("すもももももももものうち。", 1)
unlist(RMeCabC("すもももももももものうち。", 1))
#-15.6-------------------------------------------------------------------------#
# Windowsの場合
data_txt <- iconv(as.matrix(data_chap15$Content),
from = "UTF-8", to = "CP932", sub="byte")
# MacOS, Google Colabの場合
data_txt <- as.matrix(data_chap15$Content)
#-15.7-------------------------------------------------------------------------#
text_out <- unlist(RMeCabC(data_txt[1], 1))
text_out
#-15.8-------------------------------------------------------------------------#
c1 <- nchar(text_out) > 1
c2 <- names(text_out) != "記号"
c3 <- names(text_out) != "助詞"
c4 <- names(text_out) != "助動詞"
#-15.9-------------------------------------------------------------------------#
c5 <- !grepl("^[0-9]{1,}$", text_out)
c6 <- !grepl("^[0-9０-９]{1,}$", text_out)
c7 <- !grepl("^[\u3040-\u309F]{1,3}$", text_out)
#-15.10------------------------------------------------------------------------#
text_out <- unlist(RMeCabC(data_txt[1], 1))
c1 <- nchar(text_out) > 1
c2 <- names(text_out) != "記号"
c3 <- names(text_out) != "助詞"
c4 <- names(text_out) != "助動詞"
c5 <- !grepl("^[0-9]{1,}$", text_out)
c6 <- !grepl("^[0-9０-９]{1,}$", text_out)
c7 <- !grepl("^[\u3040-\u309F]{1,3}$", text_out)
text_out <- text_out[c1&c2&c3&c4&c5&c6&c7]
#-15.11------------------------------------------------------------------------#
text_out
#-15.12------------------------------------------------------------------------#
corpus_data <- NULL
for (n in 1:length(data_txt)){
text_out <- unlist(RMeCabC(data_txt[n],1))
c1 <- nchar(text_out) > 1
c2 <- names(text_out) != "記号"
c3 <- names(text_out) != "助詞"
c4 <- names(text_out) != "助動詞"
c5 <- !grepl("^[0-9]{1,}$", text_out)
c6 <- !grepl("^[0-9０-９]{1,}$", text_out)
c7 <- !grepl("^[\u3040-\u309F]{1,3}$", text_out)
text_out <- text_out[c1&c2&c3&c4&c5&c6&c7]
corpus_data <- rbind(corpus_data, cbind(n, text_out))
}
#-15.13------------------------------------------------------------------------#
# 本書の出力を再現したい場合：分解済みコーパス
corpus_data <- read.xlsx("chapter_15_corpus.xlsx")[,-1]
#-15.14------------------------------------------------------------------------#
vocab_tab <- sort(table(corpus_data[,2]), decreasing = TRUE)
#-15.15------------------------------------------------------------------------#
vocab_tab[1:50]
#-15.16------------------------------------------------------------------------#
word_list <- names(table(corpus_data[,2]))
#-15.17------------------------------------------------------------------------#
dat_x <- as.numeric(factor(corpus_data[,1],
levels = 1:nrow(data_txt)))
dat_w <- as.numeric(factor(corpus_data[,2]))
#-15.18------------------------------------------------------------------------#
N <- length(dat_x)
V <- length(word_list)
D <- nrow(data_txt)
K <- 6
#-15.19------------------------------------------------------------------------#
set.seed(1)
dat_z <- ceiling(runif(N)*K)
#-15.20------------------------------------------------------------------------#
Nkv <- matrix(0, K, V)
Nk <- numeric(K)
Ndk <- matrix(0, D, K)
Nd <- table(dat_x)
for (i in 1:N){
Nkv[dat_z[i], dat_w[i]] <- Nkv[dat_z[i], dat_w[i]] + 1
Ndk[dat_x[i], dat_z[i]] <- Ndk[dat_x[i], dat_z[i]] + 1
Nk[dat_z[i]] <- Nk[dat_z[i]] + 1
}
#-15.21------------------------------------------------------------------------#
NN <- 1500
b0 <- 0.1
a0 <- K / 50
est_z <- matrix(0, NN, N)
#-15.22------------------------------------------------------------------------#
set.seed(1)
for (nn in 1:NN){
for (i in 1:N){
Nkv[dat_z[i], dat_w[i]] <- Nkv[dat_z[i], dat_w[i]] - 1
Ndk[dat_x[i], dat_z[i]] <- Ndk[dat_x[i], dat_z[i]] - 1
Nk[dat_z[i]] <- Nk[dat_z[i]] - 1
eta <- ((Nkv[,dat_w[i]] + b0) / (Nk + V * b0)) *
(Ndk[dat_x[i],] + a0)
p1 <- eta / sum(eta)
dat_z[i] <- which(rmultinom(1, 1, p1) == 1)
Nkv[dat_z[i],dat_w[i]] <- Nkv[dat_z[i],dat_w[i]] + 1
Ndk[dat_x[i],dat_z[i]] <- Ndk[dat_x[i],dat_z[i]] + 1
Nk[dat_z[i]] <- Nk[dat_z[i]] + 1
}
est_z[nn,] <- dat_z
#print(nn)
}
#-15.23------------------------------------------------------------------------#
res_z <- dat_z
for (i in 1:N){
res_z[i] <- which.max(table(factor(est_z[501:NN,i], levels = 1:K)))
}
#-15.24------------------------------------------------------------------------#
Nkv <- matrix(0, K, V)
Nk <- numeric(K)
Ndk <- matrix(0, D, K)
for (i in 1:N){
Nkv[res_z[i], dat_w[i]] <- Nkv[res_z[i], dat_w[i]] + 1
Ndk[dat_x[i], res_z[i]] <- Ndk[dat_x[i], res_z[i]] + 1
Nk[res_z[i]] <- Nk[res_z[i]] + 1
}
#-15.25------------------------------------------------------------------------#
res_phi <- matrix(0, K, V)
for (k in 1:K){
res_phi[k,] <- (Nkv[k,] + b0) / (Nk[k] + V * b0)
}
colnames(res_phi) <- word_list
res_theta <- matrix(0, D, K)
for (d in 1:D){
res_theta[d,] <- (Ndk[d,] + a0) / (Nd[d] + K * a0)
}
#-15.26------------------------------------------------------------------------#
top_words <- matrix(0, K, 10)
for (k in 1:K){
top_words[k,] <- names(sort(res_phi[k,], decreasing = TRUE))[1:10]
}
#-15.27------------------------------------------------------------------------#
t(top_words)
#-15.28------------------------------------------------------------------------#
table(data_chap15$Satisfaction)
#-15.29------------------------------------------------------------------------#
table(data_chap15$OS)
#-15.30------------------------------------------------------------------------#
os_list <- names(sort(table(data_chap15$OS)))
os_theta <- matrix(0, length(os_list), K)
rownames(os_theta) <- os_list
for (j in 1:length(os_list)){
os_theta[j,] <- colMeans(res_theta[data_chap15$OS == os_list[j],])
}
#-15.31------------------------------------------------------------------------#
barplot(t(os_theta), beside = TRUE)
#-15.32------------------------------------------------------------------------#
year <- floor(data_chap15$Month / 12)
table(year)
year[year > 4] <- 4
#-15.33------------------------------------------------------------------------#
year_list <- names(table(year))
year_theta <- matrix(0,length(year_list), K)
rownames(year_theta) <- year_list
for (j in 1:length(year_list)){
year_theta[j,] <- colMeans(res_theta[year == year_list[j],])
}
#-15.34------------------------------------------------------------------------#
matplot(year_theta, type = "b",
pch = as.character(c(1:9, 0)),
col = "black", lty = 1, xaxt = "n", xlab = "使用年数")
axis(1, at = 1:length(year_list), labels = year_list)
#------------------------------------------------------------------------------#