PI データ分析ワークフロー
データ分析チーム
2026-04-23
はじめに
このレポートでは、PI Web APIを使用してPIシステムからデータを取得し、Rを用いて可視化する一連の流れを示します。 今回は、以下のタグのデータを取得対象とします： M2TI4405.PV, M2PI4407.PV, M2PI4410W.PV, M2PI4411W.PV, M2FC4403.PV, M2CI2901.CPV
ステップ 1: ライブラリの読み込み
APIリクエストには httr、JSONの解析には
jsonlite、データの整形・可視化には dplyr,
lubridate, ggplot2 を使用します。
library(httr)
library(jsonlite)
library(ggplot2)
library(dplyr)
library(lubridate)
library(tidyr)
library(DT)ステップ 2: 設定
PI Web APIへの接続パラメータを設定します。 ※
URL等の設定値はリポジトリ内の設定ファイル
(expressions.tmdl) を参照しています。
# PIシステムの設定
PI_WEB_API_BASE_URL <- "https://hinuta105.ccit.ad.sharedom.net/piwebapi"
PI_SERVER_NAME <- "HINUTA101C"
TAG_NAMES <- c("M2TI4405.PV", "M2PI4407.PV", "M2PI4410W.PV", "M2PI4411W.PV", "M2FC4403.PV", "M2CI2901.CPV")
# 認証設定
# 'auth_config.R' から認証情報を読み込みます
if (file.exists("auth_config.R")) {
source("auth_config.R")
auth_config <- authenticate(AUTH_USER, AUTH_PASS)
} else {
stop("認証設定ファイル 'auth_config.R' が見つかりません。'auth_config_template.R' をコピーして設定してください。")
}
ssl_config <- config(ssl_verifypeer = 0) # 社内サーバー等のためSSL検証をスキップステップ 3: 関数定義
APIと対話するためのヘルパー関数を定義します。 1.
get_web_id: タグ名からWebIDを取得する関数 2.
get_interpolated_data:
WebIDを使用して指定期間・周期の補間データ（Interpolated
Data）を取得する関数
get_web_id <- function(tag_name) {
path <- sprintf("\\\\%s\\%s", PI_SERVER_NAME, tag_name)
url <- paste0(PI_WEB_API_BASE_URL, "/points")
response <- GET(url, query = list(path = path), auth_config, ssl_config)
if (status_code(response) != 200) {
warning(paste("WebIDの取得に失敗しました:", tag_name))
return(NULL)
}
content_json <- fromJSON(content(response, "text", encoding = "UTF-8"))
return(content_json$WebId)
}
get_interpolated_data <- function(web_id, tag_name, start_time = "2021-01-01T00:00:00+09:00", end_time = "*", interval = "1d") {
url <- paste0(PI_WEB_API_BASE_URL, "/streams/", web_id, "/interpolated")
response <- GET(url, query = list(startTime = start_time, endTime = end_time, interval = interval),
auth_config, ssl_config)
if (status_code(response) != 200) {
warning(paste("データの取得に失敗しました:", tag_name))
return(NULL)
}
content_json <- fromJSON(content(response, "text", encoding = "UTF-8"))
items <- content_json$Items
if (is.null(items) || length(items) == 0) return(NULL)
# Handle complex value (e.g. Digital State object or System Error)
vals <- items$Value
if (is.data.frame(vals)) {
# If it's a data frame, we want the "Value" column if it exists, otherwise the first column.
# We must NOT unlist the entire data frame as that concatenates columns.
if ("Value" %in% names(vals)) {
vals <- vals$Value
} else {
vals <- vals[[1]]
}
} else if (is.list(vals)) {
# If it's a list of objects (e.g. from JSON), extract the Value field or the first element
# Use sapply to ensure we get a vector of the same length, handling NULLs
vals <- sapply(vals, function(x) {
if (is.list(x) || is.data.frame(x)) {
if ("Value" %in% names(x)) return(x$Value)
if (length(x) > 0) return(x[[1]])
return(NA)
} else {
return(x)
}
})
}
data.frame(
Timestamp = items$Timestamp,
Value = vals,
Tag = tag_name,
stringsAsFactors = FALSE
)
}ステップ 4: データ取得の実行
定義したタグのリストをループ処理し、WebIDを取得した後、2021年以降の日次データ（interval = "1d"）を取得します。
all_data <- data.frame()
for (tag in TAG_NAMES) {
web_id <- get_web_id(tag)
if (!is.null(web_id)) {
tag_data <- get_interpolated_data(web_id, tag)
if (!is.null(tag_data)) {
all_data <- rbind(all_data, tag_data)
}
}
}
# データの先頭を表示
head(all_data)## Timestamp Value Tag
## 1 2020-12-31T15:00:00Z 93.50242 M2TI4405.PV
## 2 2021-01-01T15:00:00Z 93.54465 M2TI4405.PV
## 3 2021-01-02T15:00:00Z 93.48495 M2TI4405.PV
## 4 2021-01-03T15:00:00Z 93.54973 M2TI4405.PV
## 5 2021-01-04T15:00:00Z 93.57597 M2TI4405.PV
## 6 2021-01-05T15:00:00Z 93.48778 M2TI4405.PVステップ 5: データ整形と可視化
取得したデータのタイムスタンプを日付型に変換し、値を数値型に変換した後、グラフを描画します。
if (nrow(all_data) > 0) {
# 型変換とタイムゾーン変換（日本時間）
all_data$Timestamp <- with_tz(ymd_hms(all_data$Timestamp), tzone = "Asia/Tokyo")
all_data$Value <- as.numeric(all_data$Value)
# 数値変換でNAになったデータ（システムステート等）を除去
all_data <- all_data %>% filter(!is.na(Value))
# プロット
ggplot(all_data, aes(x = Timestamp, y = Value, color = Tag)) +
geom_line() +
labs(title = "PIデータトレンド",
subtitle = paste(TAG_NAMES, collapse = ", "),
x = "日時",
y = "値") +
theme_minimal()
} else {
print("表示可能なデータがありません。")
}
ステップ 6: 総括伝熱係数 (U値) の算出
取得したデータから、U_Calculation_Method.md に基づいて総括伝熱係数
U を算出します。
if (nrow(all_data) > 0) {
# データをタグごとに横持ち（ワイドフォーマット）に変換
wide_data <- all_data %>%
select(Timestamp, Tag, Value) %>%
# 重複がある場合は平均をとるか最初の値をとる（日次データなので基本は1日1行）
group_by(Timestamp, Tag) %>%
summarise(Value = mean(Value, na.rm = TRUE), .groups = 'drop') %>%
pivot_wider(names_from = Tag, values_from = Value)
# 必要なタグが全て存在するか確認
required_tags <- c("M2FC4403.PV", "M2TI4405.PV", "M2PI4410W.PV", "M2PI4407.PV", "M2PI4411W.PV")
if (all(required_tags %in% names(wide_data))) {
# U値の計算プロセス
calc_data <- wide_data %>%
mutate(
# ステップ 1: 熱量 q の算出
q = M2FC4403.PV * 2257.05 * 1000 / 3600 / 1000,
# ステップ 2: ポンプ流量 pump_flow の算出
pump_inlet = (M2PI4407.PV / 760 * 101.325) + (950 * 9.81 * (1.4 + 6.125) / 1000),
pump_outlet = M2PI4411W.PV + 101.325,
lifting_height = (pump_outlet - pump_inlet) * 1000 / 9.81 / 950,
# 二次方程式の適用
pump_flow = (-0.6857 * lifting_height^2 + 0.8571 * lifting_height + 43.589) * 60 * 950 / 1000,
# ステップ 3: 温度差の算出
t1 = M2TI4405.PV,
t2 = t1 + q * 1000 / (pump_flow * 1000 / 3600) / 2.09,
# LN は R では log()
t3_t4 = 27.838 * log(M2PI4410W.PV + 101.325) - 28.559,
# ステップ 4: 対数平均温度差 delta_tlm (LMTD)
dT1 = t3_t4 - t1,
dT2 = t3_t4 - t2,
delta_tlm = ifelse(
dT1 > 0 & dT2 > 0 & dT1 != dT2,
(dT1 - dT2) / log(dT1 / dT2),
NA # 計算不可またはゼロ割りの場合
),
# 最終的な U の算出
U = ifelse(
!is.na(delta_tlm) & delta_tlm != 0,
(q * 1000 * 1000 / 367) / delta_tlm,
NA
)
) %>%
select(Timestamp, M2PI4410W = M2PI4410W.PV, q, pump_flow, t1, t2, t3_t4, delta_tlm, U) %>%
arrange(desc(Timestamp)) # 最新のデータを上に
# 計算結果の一部をコンソールに表示
head(calc_data)
} else {
print("U値の計算に必要なタグのデータが揃っていません。")
}
}## # A tibble: 6 × 9
## Timestamp M2PI4410W q pump_flow t1 t2 t3_t4 delta_tlm
## <dttm> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2026-04-23 00:00:00 -0.771 0 2388. 18.7 18.7 99.8 NA
## 2 2026-04-22 00:00:00 -0.781 0 1127. 24.4 24.4 99.8 NA
## 3 2026-04-21 00:00:00 -0.678 0 2497. 25.3 25.3 99.8 NA
## 4 2026-04-20 00:00:00 -0.701 0 1014. 35.5 35.5 99.8 NA
## 5 2026-04-19 00:00:00 -0.699 0 1017. 39.7 39.7 99.8 NA
## 6 2026-04-18 00:00:00 -0.715 0 1023. 45.5 45.5 99.8 NA
## # ℹ 1 more variable: U <dbl>ステップ 7: U値のテーブル表示
計算された U値 および各種中間指標を DT::datatable
を用いて見やすいテーブルとして表示します。
if (exists("calc_data") && nrow(calc_data) > 0) {
datatable(
calc_data,
options = list(
pageLength = 10,
autoWidth = TRUE,
scrollX = TRUE,
dom = 'Bfrtip'
),
caption = "総括伝熱係数 (U) および関連パラメータの計算結果",
rownames = FALSE
) %>%
formatRound(columns = c("M2PI4410W", "q", "pump_flow", "t1", "t2", "t3_t4", "delta_tlm", "U"), digits = 2) %>%
formatDate("Timestamp", method = "toLocaleString")
} else {
print("表示できる計算結果がありません。")
}ステップ 8: キャンペーン期間ごとの最大値取得と詰まり本数の評価
Campaign_Evaluation_Methodology.md
に基づき、各キャンペーンの「初期」と「末期」における最大ジャケット圧力時のデータを抽出し、U値との相関を確認します。
if (exists("calc_data") && nrow(calc_data) > 0 && file.exists("campaign_periods.csv")) {
# キャンペーン期間データの読み込み
campaigns <- read.csv("campaign_periods.csv", stringsAsFactors = FALSE)
campaigns$Start_Date <- as.Date(campaigns$Start_Date)
campaigns$End_Date <- as.Date(campaigns$End_Date)
campaign_results <- data.frame()
for (i in 1:nrow(campaigns)) {
period <- campaigns$Period[i]
start_d <- campaigns$Start_Date[i]
end_d <- campaigns$End_Date[i]
if (is.na(start_d)) next
# 初期 (Beginning): 開始日から1ヶ月間
beg_end_d <- start_d + months(1)
beg_data <- calc_data %>% filter(as.Date(Timestamp) >= start_d & as.Date(Timestamp) <= beg_end_d)
if (nrow(beg_data) > 0) {
max_row <- beg_data[which.max(beg_data$M2PI4410W), ]
# Use check.names=FALSE style or matching the literal name
col_beg <- grep("詰まり本数_前半", names(campaigns), value=TRUE)[1]
nos_val <- as.numeric(campaigns[i, col_beg])
campaign_results <- rbind(campaign_results, data.frame(
Period = period,
Type = "Beginning",
PeriodType = paste0(period, "_Beginning"),
Timestamp = max_row$Timestamp,
Max_Jacket_P = max_row$M2PI4410W,
U_Value = max_row$U,
Flow_Rate = max_row$pump_flow,
NoS = nos_val
))
}
# 末期 (End): 終了日の1ヶ月前から終了日まで
if (!is.na(end_d)) {
end_start_d <- end_d - months(1)
end_data <- calc_data %>% filter(as.Date(Timestamp) >= end_start_d & as.Date(Timestamp) <= end_d)
if (nrow(end_data) > 0) {
max_row <- end_data[which.max(end_data$M2PI4410W), ]
col_end <- grep("詰まり本数_後半", names(campaigns), value=TRUE)[1]
nos_val <- as.numeric(campaigns[i, col_end])
campaign_results <- rbind(campaign_results, data.frame(
Period = period,
Type = "End",
PeriodType = paste0(period, "_End"),
Timestamp = max_row$Timestamp,
Max_Jacket_P = max_row$M2PI4410W,
U_Value = max_row$U,
Flow_Rate = max_row$pump_flow,
NoS = nos_val
))
}
}
}
# 結果のテーブル表示
if (nrow(campaign_results) > 0) {
library(htmltools)
dt <- datatable(
campaign_results,
options = list(pageLength = 10, autoWidth = TRUE, scrollX = TRUE),
caption = "キャンペーン期間ごとの最大ジャケット圧力抽出結果",
rownames = FALSE
) %>%
formatRound(columns = c("Max_Jacket_P", "U_Value"), digits = 2) %>%
formatDate("Timestamp", method = "toLocaleString")
print(tagList(dt))
# 散布図の作成 1 (最大ジャケット圧力 vs U値)
p_campaign_u <- ggplot(campaign_results, aes(x = Max_Jacket_P, y = U_Value)) +
geom_point(aes(color = Type, size = NoS), alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE, color = "blue", fill = "lightblue", alpha = 0.3) +
labs(
title = "キャンペーン期間ごとの最大ジャケット圧力とU値の相関",
subtitle = "各キャンペーンの初期(1ヶ月)・末期(1ヶ月)における最大ジャケット圧力時のデータ",
x = "最大ジャケット圧力 (M2PI4410W.PV)",
y = "総括伝熱係数 (U値)",
color = "期間",
size = "詰まり本数(NoS)"
) +
theme_minimal()
print(p_campaign_u)
# 相関分析結果を保存するリスト
model_stats <- data.frame(
Parameter = character(),
R2 = numeric(),
p_value = numeric(),
Equation = character(),
stringsAsFactors = FALSE
)
# 散布図の作成 2 (最大ジャケット圧力 vs 詰まり本数)
valid_nos_p <- campaign_results %>% filter(!is.na(NoS), !is.na(Max_Jacket_P))
if(nrow(valid_nos_p) > 2) {
model_p <- lm(NoS ~ Max_Jacket_P, data = valid_nos_p)
sum_p <- summary(model_p)
r2_p <- sum_p$r.squared
pval_p <- sum_p$coefficients[2, 4]
eq_p <- sprintf("y = %.2fx %s %.2f", coef(model_p)[2], ifelse(coef(model_p)[1] >= 0, "+", "-"), abs(coef(model_p)[1]))
stats_str_p <- sprintf("回帰式: %s | R^2 = %.3f | p-value = %.3f", eq_p, r2_p, pval_p)
model_stats <- rbind(model_stats, data.frame(Parameter = "ジャケット圧力", R2 = r2_p, p_value = pval_p, Equation = eq_p))
p_campaign_nos_p <- ggplot(valid_nos_p, aes(x = Max_Jacket_P, y = NoS)) +
geom_point(aes(color = Type), size = 3, alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE, color = "darkred", fill = "pink", alpha = 0.3) +
labs(
title = "キャンペーン期間ごとの最大ジャケット圧力と詰まり本数の相関",
subtitle = paste("各キャンペーンの初期(1ヶ月)・末期(1ヶ月)における最大ジャケット圧力時のデータ\n", stats_str_p),
x = "最大ジャケット圧力 (M2PI4410W.PV)",
y = "詰まり本数 (NoS)",
color = "期間"
) +
theme_minimal()
print(p_campaign_nos_p)
}
# 散布図の作成 3 (U値 vs 詰まり本数)
valid_nos_u <- campaign_results %>% filter(!is.na(NoS), !is.na(U_Value))
if(nrow(valid_nos_u) > 2) {
model_u <- lm(NoS ~ U_Value, data = valid_nos_u)
sum_u <- summary(model_u)
r2_u <- sum_u$r.squared
pval_u <- sum_u$coefficients[2, 4]
eq_u <- sprintf("y = %.2fx %s %.2f", coef(model_u)[2], ifelse(coef(model_u)[1] >= 0, "+", "-"), abs(coef(model_u)[1]))
stats_str_u <- sprintf("回帰式: %s | R^2 = %.3f | p-value = %.3f", eq_u, r2_u, pval_u)
model_stats <- rbind(model_stats, data.frame(Parameter = "U値", R2 = r2_u, p_value = pval_u, Equation = eq_u))
p_campaign_nos_u <- ggplot(valid_nos_u, aes(x = U_Value, y = NoS)) +
geom_point(aes(color = Type), size = 3, alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE, color = "darkgreen", fill = "lightgreen", alpha = 0.3) +
labs(
title = "キャンペーン期間ごとのU値と詰まり本数の相関",
subtitle = paste("各キャンペーンの初期(1ヶ月)・末期(1ヶ月)における最大ジャケット圧力時のデータ\n", stats_str_u),
x = "総括伝熱係数 (U値)",
y = "詰まり本数 (NoS)",
color = "期間"
) +
theme_minimal()
print(p_campaign_nos_u)
}
# 散布図の作成 4 (ポンプ流量 vs 詰まり本数)
valid_nos_f <- campaign_results %>% filter(!is.na(NoS), !is.na(Flow_Rate))
if(nrow(valid_nos_f) > 2) {
model_f <- lm(NoS ~ Flow_Rate, data = valid_nos_f)
sum_f <- summary(model_f)
r2_f <- sum_f$r.squared
pval_f <- sum_f$coefficients[2, 4]
eq_f <- sprintf("y = %.2fx %s %.2f", coef(model_f)[2], ifelse(coef(model_f)[1] >= 0, "+", "-"), abs(coef(model_f)[1]))
stats_str_f <- sprintf("回帰式: %s | R^2 = %.3f | p-value = %.3f", eq_f, r2_f, pval_f)
model_stats <- rbind(model_stats, data.frame(Parameter = "ポンプ流量", R2 = r2_f, p_value = pval_f, Equation = eq_f))
p_campaign_nos_f <- ggplot(valid_nos_f, aes(x = Flow_Rate, y = NoS)) +
geom_point(aes(color = Type), size = 3, alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE, color = "purple", fill = "thistle", alpha = 0.3) +
labs(
title = "キャンペーン期間ごとのポンプ流量と詰まり本数の相関",
subtitle = paste("各キャンペーンの初期(1ヶ月)・末期(1ヶ月)における最大ジャケット圧力時のデータ\n", stats_str_f),
x = "ポンプ流量 (pump_flow)",
y = "詰まり本数 (NoS)",
color = "期間"
) +
theme_minimal()
print(p_campaign_nos_f)
}
# 詰まり本数との相関モデル比較の表示
if (nrow(model_stats) > 0) {
model_stats <- model_stats %>% arrange(desc(R2))
print(datatable(
model_stats,
options = list(pageLength = 5, autoWidth = TRUE, dom = 't'),
caption = "詰まり本数との相関強度の比較 (R^2 降順)",
rownames = FALSE
) %>%
formatRound(columns = c("R2", "p_value"), digits = 3))
}
} else {
print("キャンペーン評価用のデータが抽出できませんでした。")
}
} else {
print("キャンペーン評価に必要なデータ(calc_data)または campaign_periods.csv が存在しません。")
}## <div class="datatables html-widget html-fill-item" id="htmlwidget-b8384fd0f1eb692f436c" style="width:100%;height:auto;"></div>
## <script type="application/json" data-for="htmlwidget-b8384fd0f1eb692f436c">{"x":{"filter":"none","vertical":false,"caption":"<caption>キャンペーン期間ごとの最大ジャケット圧力抽出結果<\/caption>","data":[[2103,2103,2109,2109,2206,2206,2302,2302,2310,2310,2404,2404,2502,2502],["Beginning","End","Beginning","End","Beginning","End","Beginning","End","Beginning","End","Beginning","End","Beginning","End"],["2103_Beginning","2103_End","2109_Beginning","2109_End","2206_Beginning","2206_End","2302_Beginning","2302_End","2310_Beginning","2310_End","2404_Beginning","2404_End","2502_Beginning","2502_End"],["2021-04-13T15:00:00Z","2021-09-08T15:00:00Z","2021-10-26T15:00:00Z","2022-02-06T15:00:00Z","2022-07-26T15:00:00Z","2023-01-13T15:00:00Z","2023-03-15T15:00:00Z","2023-09-18T15:00:00Z","2023-11-17T15:00:00Z","2024-03-12T15:00:00Z","2024-04-14T15:00:00Z","2025-06-27T15:00:00Z","2025-07-30T15:00:00Z","2025-12-23T15:00:00Z"],[12.09471,13.94503,4.7961,18.06892,3.6427,-0.5194795,-1.391797,-0.5811917,-0.6167568,-0.9513991000000001,-1.175398,3.71782637,1.37394071,1.12084067],[448.6526527294358,372.7256412851982,526.9320277444729,351.2972468028649,433.386667005834,519.1924538640029,null,464.490907696674,null,null,null,480.9218610433566,479.5211867889088,593.7199198328791],[1925.796003595965,1570.825104049807,1804.486510953861,1771.142357478077,1780.787263384056,2369.987806467864,2079.520604753937,2245.628717555061,2097.86518490855,98.31618093169857,2435.517765740384,2385.43878695637,2139.425577945138,2009.563149209568],[138,196,123,162,0,9,9,18,18,26,null,null,26,null]],"container":"<table class=\"display\">\n <thead>\n <tr>\n <th>Period<\/th>\n <th>Type<\/th>\n <th>PeriodType<\/th>\n <th>Timestamp<\/th>\n <th>Max_Jacket_P<\/th>\n <th>U_Value<\/th>\n <th>Flow_Rate<\/th>\n <th>NoS<\/th>\n <\/tr>\n <\/thead>\n<\/table>","options":{"pageLength":10,"autoWidth":true,"scrollX":true,"columnDefs":[{"targets":3,"render":"function(data, type, row, meta) {\n return type !== 'display' ? data : DTWidget.formatDate(data, \"toLocaleString\");\n }"},{"targets":4,"render":"function(data, type, row, meta) {\n return type !== 'display' ? data : DTWidget.formatRound(data, 2, 3, \",\", \".\", null);\n }"},{"targets":5,"render":"function(data, type, row, meta) {\n return type !== 'display' ? data : DTWidget.formatRound(data, 2, 3, \",\", \".\", null);\n }"},{"className":"dt-right","targets":[0,4,5,6,7]},{"name":"Period","targets":0},{"name":"Type","targets":1},{"name":"PeriodType","targets":2},{"name":"Timestamp","targets":3},{"name":"Max_Jacket_P","targets":4},{"name":"U_Value","targets":5},{"name":"Flow_Rate","targets":6},{"name":"NoS","targets":7}],"order":[],"orderClasses":false}},"evals":["options.columnDefs.0.render","options.columnDefs.1.render","options.columnDefs.2.render"],"jsHooks":[]}</script>
ステップ 9: 予測モデルの構築と評価 (2502期の予測)
2310以前のデータを訓練データ（学習データ）として線形回帰モデルを構築し、2502期（テストデータ）の詰まり本数を予測します。実際の2502期の詰まり本数（48本）と比較し、モデルの予測精度を評価します。
if (exists("campaign_results") && nrow(campaign_results) > 0) {
# 訓練データ: Periodが2310以下、かつNoSに欠損値がない行
train_data <- campaign_results %>%
filter(as.numeric(Period) <= 2310, !is.na(NoS))
# テストデータの抽出: 2502期のデータ
test_data_base <- campaign_results %>% filter(Period == "2502")
# テストデータの予測変数を取得 (初期と末期で分かれている場合は、予測したいフェーズに合わせますが、ここでは最大値を持つ行を採用するか、テストデータ群として扱います。依頼の意図としては2502期のある特定時点での予測値と実測値の比較と考えられますので、テストデータとして利用可能な行ごとに予測を行います。)
if (nrow(train_data) > 2 && nrow(test_data_base) > 0) {
# 実際の値
actual_nos_2502 <- 47
# 訓練データでのモデル構築
mod_p <- lm(NoS ~ Max_Jacket_P, data = train_data)
mod_u <- lm(NoS ~ U_Value, data = train_data)
mod_f <- lm(NoS ~ Flow_Rate, data = train_data)
eval_results <- data.frame()
for(i in 1:nrow(test_data_base)) {
row_data <- test_data_base[i, ]
# 予測の実行
pred_p <- predict(mod_p, newdata = row_data)
pred_u <- predict(mod_u, newdata = row_data)
pred_f <- predict(mod_f, newdata = row_data)
# 評価指標 (絶対誤差と絶対誤差率 MAPE)
error_p <- abs(pred_p - actual_nos_2502)
error_u <- abs(pred_u - actual_nos_2502)
error_f <- abs(pred_f - actual_nos_2502)
mape_p <- (error_p / actual_nos_2502) * 100
mape_u <- (error_u / actual_nos_2502) * 100
mape_f <- (error_f / actual_nos_2502) * 100
eval_results <- rbind(eval_results, data.frame(
PeriodType = row_data$PeriodType,
Actual_NoS = actual_nos_2502,
Predictor = c("ジャケット圧力", "U値", "ポンプ流量"),
Predicted_NoS = c(pred_p, pred_u, pred_f),
Absolute_Error = c(error_p, error_u, error_f),
MAPE_Percent = c(mape_p, mape_u, mape_f)
))
}
# 評価結果の列を追加
eval_results <- eval_results %>% mutate(
Evaluation = case_when(
MAPE_Percent < 10 ~ "非常に良い予測",
MAPE_Percent >= 10 & MAPE_Percent < 20 ~ "良い予測",
MAPE_Percent >= 20 & MAPE_Percent < 50 ~ "妥当な予測",
TRUE ~ "不正確な予測"
)
)
# 列名を日本語に変更してわかりやすくする
eval_results_display <- eval_results %>%
rename(
`期間` = PeriodType,
`実際の詰まり本数` = Actual_NoS,
`予測に使用した指標` = Predictor,
`予想本数` = Predicted_NoS,
`誤差（絶対値）` = Absolute_Error,
`誤差率 (MAPE %)` = MAPE_Percent,
`評価結果` = Evaluation
)
# MAPEで昇順ソートして最も精度が良いものを上に
eval_results_display <- eval_results_display %>% arrange(`誤差率 (MAPE %)`)
# 評価コメントの自動生成 (表の前に出力)
best_model <- eval_results %>% arrange(MAPE_Percent) %>% slice(1)
cat(sprintf("\n**評価コメント:**\n一般的な評価指標（MAPE: 平均絶対パーセント誤差）に照らすと、%sの時の '%s' を使ったモデルが最も誤差が少なく、予測値は約 %.1f 本（誤差 %.1f%%、評価: %s）となりました。\n※ 評価基準: MAPEが10%%未満で「非常に良い」、10%%〜20%%未満で「良い」、20%%〜50%%未満で「妥当」、50%%以上で「不正確」としています。\n\n",
best_model$PeriodType, best_model$Predictor, best_model$Predicted_NoS, best_model$MAPE_Percent, best_model$Evaluation))
datatable(
eval_results_display,
options = list(pageLength = 10, autoWidth = TRUE, scrollX = TRUE),
caption = "2502期の予測モデル評価結果 (実際の詰まり本数: 47本)",
rownames = FALSE
) %>%
formatRound(columns = c("予想本数", "誤差（絶対値）", "誤差率 (MAPE %)"), digits = 2)
} else {
print("モデルの訓練に必要な2310以前のデータ、または2502期のテストデータが不足しています。")
}
}##
## **評価コメント:**
## 一般的な評価指標（MAPE: 平均絶対パーセント誤差）に照らすと、2502_Beginningの時の 'ジャケット圧力' を使ったモデルが最も誤差が少なく、予測値は約 37.0 本（誤差 21.2%、評価: 妥当な予測）となりました。
## ※ 評価基準: MAPEが10%未満で「非常に良い」、10%〜20%未満で「良い」、20%〜50%未満で「妥当」、50%以上で「不正確」としています。