#VGAMのインストールが済んでいない場合，以下の行を実行してください．

install.packages("VGAM")

#パッケージの読み込み

library(VGAM)

#以下を実行したときに表示されるフォルダ名にデータ（v08c5tv.csv）を置いてください．

getwd()

#データの読み込み

c5tv <- read.csv("v08c5tv.csv")

#問2

#1変数多項ロジットモデル

c5logitm1<-vglm(購入機種 ~ 基本性能重視度, data = c5tv, family = multinomial())

#カテゴリの並び順を確認(最終カテゴリが魅力を1と調整した基準．この場合はEが基準，1がA，2がB…)

head(c5logitm1@y)

#基本性能重視度4〜20のデータフレーム作成．

c5newdata1 = data.frame("基本性能重視度" = c(4:20))

#4〜20それぞれの場合の予測値の算出

c5predict1 <- predict(c5logitm1, c5newdata1, type = "response")

#予測値の行名を4〜20にする

rownames(c5predict1) <- t(c5newdata1["基本性能重視度"])

#グラフ描画

#まず空のグラフでグラフの外枠だけ作る．

plot(0, 0, type = "n", xlim = c(4, 20), ylim = c(0, 0.6),xlab = "基本性能重視度", ylab = "購入確率")

#A〜Eのグラフの色設定

cols <- c("black", "red", "blue", "orange", "green")

#1列目から5列目までを曲線として書き加える．

for (i in 1:ncol(c5predict1)){

lines(rownames(c5predict1), c5predict1[, i], col = cols[i])

}

#グラフに凡例を加える

legend("top", legend = colnames(c5predict1), col = cols, lty = 1)

#問3

#10,11の確率

predict(c5logitm1, data.frame("基本性能重視度" = c(10, 11)), type = "response")

#10,11のオッズ．log(mu[,1]/mu[,5]) の1が基本性能重視度10の場合のAのオッズ

(c5odds <- exp(predict(c5logitm1, data.frame("基本性能重視度" = c(10, 11)))))

#オッズ比

c5odds[2,1]/c5odds[1,1]

#重みを使ったオッズ比計算．coef(分析結果)で係数を取り出すことができる．この部分はglmと使い方が異なる．

#すべての係数が並んだベクトルとして取り出され，その5番目がAの重み．exp(重み)オッズ比．

exp(coef(c5logitm1)[5])

#問4

#2変数多項ロジットモデル．共用が1の場合共用，0の場合個人用．

c5logitm2<-vglm(購入機種 ~ 基本性能重視度 + 共用, data = c5tv, family = multinomial())

#問5

#共用の場合と個人用の場合，それぞにおける基本性能重視度10の時の購入確率．

predict(c5logitm2, data.frame("共用" = c(0, 1),"基本性能重視度" = c(10, 10)), type = "response")

coef(c5logitm1)

#====================================================

#2022.09.09 齋藤朗宏・荘島宏二郎