# 必要なライブラリのインストール（初回のみ）
# install.packages("ggplot2")
# install.packages("lmtest") # 尤度比検定（lrtest）などに使用

# ライブラリのロード
library(ggplot2) # グラフ描画
library(lmtest)  # 尤度比検定（オプション）

# --- 準備：InputPathの設定 ---
# ここにご自身のCSVファイルのフルパスを指定してください
InputPath <- "C:/data/my_logistic_data.csv" # 例：Windows
# InputPath <- "/Users/yourname/data/my_data.csv" # 例：Mac/Linux

# -----------------------------------------------------------------
# ① CSVファイルの読み込み
# -----------------------------------------------------------------

# InputPathにファイルが存在するかチェック
if (file.exists(InputPath)) {
  data <- read.csv(InputPath)
  print(paste(InputPath, "を読み込みました。"))
} else {
  # ファイルが存在しない場合、警告を出し、サンプルデータを作成
  warning(paste("警告:", InputPath, "が見つかりません。代わりにサンプルデータを生成します。"))
  
  set.seed(42) # 結果の再現性のため
  N <- 200     # データ数
  x <- rnorm(N, 5, 2) # 独立変数 x (平均5, 標準偏差2の正規分布)
  
  # ロジット（対数オッズ）を定義
  # 真のモデル: logit(p) = -5 + 1.2*x
  logit_p <- -5 + 1.2 * x
  
  # 確率 p に変換
  p <- 1 / (1 + exp(-logit_p))
  
  # 0か1の従属変数 y を生成
  y <- rbinom(N, 1, p)
  
  data <- data.frame(x = x, y = y)
  
  print("サンプルデータを生成しました。")
}

# データの先頭部分を表示
# print(head(data))

# -----------------------------------------------------------------
# ② ロジスティック回帰の実行
# -----------------------------------------------------------------
# glm (Generalized Linear Model) 関数を使用
# family = binomial(link = "logit") でロジスティック回帰を指定

# 予測変数(x)を含むフルモデル
model_full <- glm(y ~ x, data = data, family = binomial(link = "logit"))

# 比較用：切片のみの帰無モデル（Nullモデル）
model_null <- glm(y ~ 1, data = data, family = binomial(link = "logit"))

# -----------------------------------------------------------------
# ③ 結果の表示 と ④ Devianceの表示
# -----------------------------------------------------------------
# summary() を使うと、結果、ワルド検定、Devianceが一度に表示されます

summary_model <- summary(model_full)

print("--- ロジスティック回帰モデルの結果概要 ---")
print(summary_model)

cat("\n--- Deviance情報 (概要から抜粋) ---\n")
# Null Deviance: 帰無モデル（切片のみ）の逸脱度
cat("Null Deviance (帰無モデル): ", summary_model$null.deviance, 
    " (自由度: ", summary_model$df.null, ")\n")

# Residual Deviance: 予測変数(x)を含むモデルの逸脱度
cat("Residual Deviance (残差): ", summary_model$deviance, 
    " (自由度: ", summary_model$df.residual, ")\n")
cat("AIC (赤池情報量基準): ", summary_model$aic, "\n")


# -----------------------------------------------------------------
# ③-2 グラフ化 (ggplot2を使用)
# -----------------------------------------------------------------
print("--- 予測確率のグラフを生成 ---")

# 予測確率 P(y=1 | x) をデータフレームに追加
data$predicted_prob <- predict(model_full, type = "response")

# グラフを描画
plot_logistic <- ggplot(data, aes(x = x, y = y)) +
  # 観測値 (0と1) を点でプロット (jitterで点を少し散らす)
  geom_point(aes(color = as.factor(y)), 
             alpha = 0.6, 
             position = position_jitter(height = 0.03, width = 0)) +
  # ロジスティック回帰曲線（予測確率）
  geom_line(aes(y = predicted_prob), color = "blue", linewidth = 1.2) +
  labs(title = "ロジスティック回帰曲線と観測データ",
       x = "x (独立変数)",
       y = "y (観測値 / 予測確率)",
       color = "観測値 (y)") +
  theme_minimal() +
  scale_color_manual(values = c("0" = "darkorange", "1" = "steelblue")) # 色の指定

# グラフを表示
print(plot_logistic)


# -----------------------------------------------------------------
# ⑤ 尤度比検定、ワルド検定、スコア検定
# -----------------------------------------------------------------
print("\n--- 統計検定の結果 ---")

# 1. ワルド検定 (Wald Test)
# --------------------------
# ワルド検定は、`summary(model_full)` の結果に標準で含まれています。
# `Coefficients:` の `z value` と `Pr(>|z|)` が該当します。
print("--- 1. ワルド検定 (Wald Test) ---")
cat("`summary(model_full)` の係数表 (Coefficients) を参照してください。\n")
print(summary_model$coefficients)


# 2. 尤度比検定 (Likelihood Ratio Test, LRT)
# --------------------------
# 帰無モデルとフルモデルを比較します。
print("--- 2. 尤度比検定 (Likelihood Ratio Test) ---")
# statsパッケージのanova関数 ( test="LRT" を指定 )
lrt_result <- anova(model_null, model_full, test = "LRT")
print(lrt_result)
# もしくは lmtest::lrtest を使用
# print(lmtest::lrtest(model_null, model_full))


# 3. スコア検定 (Score Test or Rao Test)
# --------------------------
# フルモデルから項を削除する際の検定として `drop1` を使用します。
print("--- 3. スコア検定 (Rao Test) ---")
# stats::drop1 関数で `test = "Rao"` を指定
score_test_result <- drop1(model_full, test = "Rao")
print(score_test_result)