# try_lm02.rb.txt  ruby script (encoding: Windows-31J)
  # 単回帰分析結果のレポート(docxとして出力)
  # Web上のExcelファイルを読み込んで分析する
  # 散布図と回帰直線の図を作成
  # MarkDown原稿部分は、実は try_lm01.rb.txt と全く同じ
require "rpt"
require "erb"
require "spreadsheet"

        ## 分析の素材情報
url = "http://www.qmss.jp/databank/preliminaries/data/1-2b-2.xls" # 素材データ
sheets = Rrx::xls2ary(url)  # Excelデータを配列に変換
tcsv = Rrx::temp_make(sheets[0], ".csv")  # 第１ワークシートをcsvとして出力
ddd = tcsv.path  # csvファイル(テンポラリファイル)のフルパス名
xxx = "BMI"  # 説明変数(推定の材料)
yyy = "体脂肪率"  # 目的変数(推定したい項目)
fff_png1 = "z_lm02a.png"  # グラフ画像ファイル名その１
fff_png2 = "z_lm02b.png"  # グラフ画像ファイル名その２
File.unlink(fff_png1)  if test(?e, fff_png1)  # 画像ファイルが既にあるなら削除
File.unlink(fff_png2)  if test(?e, fff_png2)

        ## Rプログラムとmarkdown原稿の整理
src_data = DATA.read  # Rプログラムとmarkdownの雛形文書
src_data = Rrx::erb_three(src_data)
r_prog, mkd_str = src_data.split(/^###+.*?\n+/)  # Rとmarkdownの分離

        ## Rの実行
r_prog = ERB.new(r_prog).result(binding)
hs = Rrx::rexec(r_prog)  # Rの実行
tcsv.unlink  # csvテンポラリファイルの削除

        ## 分析結果の整理
cor = hs["相関"]  # 相関の検定結果
cor_result = sprintf("r=%g, df=%d, p-value=%g", cor["estimate"],
    cor["parameter"], cor["p.value"])
lm = hs["回帰"]  # 回帰分析の結果
aaa = lm["coefficients"][0]  # y切片
bbb = lm["coefficients"][1]  # xの傾き
lm2 = hs["回帰要約"]  # 回帰分析結果のsummary
rrr = lm2["r.squared"]  # 決定係数
rrr_main = sprintf("%.1f", rrr*100.0)
rrr_other = sprintf("%.1f", 100.0-rrr_main.to_f)
coef = lm2["coefficients"].dup  # 回帰係数
coef[0] = [nil, "見積り", "標準誤差", "t検定値", "t検定時のp値"]
coef[1][0] = "定数項"
        ## pandocによるdocxの生成
mkd_str = ERB.new(mkd_str).result(binding)
docx_str = Rrx::pandoc_docx(mkd_str)
File.open("z_lm02.docx", "wb") {|ff|  ff.write docx_str}
File.unlink(fff_png1)  if test(?e, fff_png1)  # 画像ファイルを削除
File.unlink(fff_png2)  if test(?e, fff_png2)

__END__
    ## Rプログラム　ここから
dtf <- read.csv("ddd", header=T)
res <- cor.test(dtf$xxx, dtf$yyy)  # 相関の検定
robj(res, "相関")
      # 散布図の作成
png(file = "fff_png1")
par(family = "Japan1GothicBBB")  # 日本語フォントにゴチックを指定
par(mar = c(3, 3, 3, 3))  # 余白
plot(dtf$xxx, dtf$yyy, type="p",
    main="図１：散布図", cex.main=2,  # メインタイトルとその大きさ
    xlab="xxx", ylab="yyy", cex.lab=2)
dev.off()
      # 回帰分析
res <- lm(yyy ~ xxx, data = dtf)
robj(res, "回帰")
res2 <- summary(res)
      # F検定のp値を算出
fx <- res2$fstatistic
res2[["f_p.value"]] <- 1 - pf(fx["value"], fx["numdf"], fx["dendf"])
robj(res2, "回帰要約")
      # 散布図と回帰直線を描く
png(file = "fff_png2")
par(family = "Japan1GothicBBB")  # 日本語フォントにゴチックを指定
par(mar = c(3, 3, 3, 3))  # 余白
plot(dtf$xxx, dtf$yyy, type="p",
    main="図２：散布図と回帰直線", cex.main=2,  # メインタイトルとその大きさ
    xlab="xxx", ylab="yyy", cex.lab=2)
abline(coef = res$coefficients, col = "blue")
dev.off()

###### 以下はmarkdownの原稿

# xxxとyyyの関係に関する回帰分析

ここにxxxとyyyを列記したデータがある。これについて回帰分析を行う。

xxxを材料にしてyyyを推定できるか調べる。
このような推定を専門的表現で「yyyをxxxに回帰させる」というらしい。

## 相関の確認

まず、回帰分析を行う意味のあるデータなのかどうかを確認する。

２つの項目の相関が認められなければ、
回帰分析の結果が実質的な意味を持たない。

「xxx」と「yyy」の相関係数を確認すると次のとおり。

<%= cor_result %>

有意な強い相関が認められる。

ちなみに、散布図は図１のようになる。

![図１：散布図](fff_png1)

----------------

## 回帰分析

強い相関が認められるので、回帰分析を行ってみた。その結果は次のとおり。

- y切片： aaa
- xの傾き： bbb

つまり、次のような推定のための計算式(回帰式)が求められる。

「yyy」＝aaa＋bbb×「xxx」

上の回帰式によって描かれる直線が回帰直線である。

散布図に回帰直線を書き入れてみると図２のようになる。

![図２：散布図と回帰直線](fff_png2)

個々の実測値と回帰式から得られる値との間には誤差(残差)がある。
すべての実測値についてこの誤差の平方和が最小になるように
考え出されたのが上の回帰式である。

----------------

## 回帰式の信頼性：決定係数

先の回帰式がどの程度信頼できるかをみる一つの手がかりは、決定係数である。

- 決定係数: rrr
- Ｆ検定の値: <%= lm2["fstatistic"][0] %> (決定係数に関するＦ検定量)
- Ｆ検定のp値: <%= lm2["f_p.value"] %>

決定係数が１に近いほど「xxx」と「yyy」の関連性が強い。

今回の rrr という値は、「yyy」の変動が「xxx」の変動によって
rrr_main％だけ説明できることを意味している。

残りのrrr_other％は、回帰式に含まれていない別の要素が関係しているとみられる。

決定係数の信頼性はＦ検定の結果により判断する。もっともっと幅広く
調査したとき、もしかすると決定係数が０になるかもしれない。その危険性を
示すのがＦ検定のp値である。これが0.05(5%)とか0.01、あるいは0.001より
小さければ、そのレベルに応じた有意性が認められることになる。

このp値が0.05以上だと、「母集団において決定係数が０である」
という仮説を、母集団に関して棄却できないことになる。

## 回帰式の信頼性：回帰係数

回帰係数(coefficient：y切片やxの傾きの関連)
に関する詳細な情報を掲げると次のとおり。

<%= Rrx::matrix2table(coef) %>

この表で、１行目の「定数項」の「見積り」はy切片の値であり、
２行目の「<%= coef[2][0] %>」の「見積り」の欄はxの傾きの値である。

重回帰分析の場合は、説明変数の個数に応じて、
３行目・４行目・……が表示されることになる。

表の右端の「t検定時のp値」は、該当の変数が実は影響力を
持たないかもしれない危険性(確率)を示す。

xの傾きが０であれば、xの値が何であってもyの値に影響しない。

xxxの傾きがbbbと見積もられてはいるが、もっともっと調査対象を拡げると、
実は０ということもあるかもしれない。その「かもしれない」の確率がp値である。
p値が十分小さければ、有意性が認められることになる。

定数項のp値については、実際に言及することは少ないと思うが、
もっともっと幅広く調査した時にy切片が０であるかもしれない確率を示す。