〜 markdownを利用してhtml, LaTeX→pdf作成 〜
最終更新日: 2013/11/20
Rとrubyの連携により、markdownの原稿を html と LaTeX の２種類に変換します。
回帰分析の簡単なレポート作成を材料として取り上げます。
回帰分析の結果について、簡単なレポートを作成します。
統計解析ソフトR、ruby、そしてrubyのライブラリkramdownを用います。拙作rrx(rubyからRを利用するためのライブラリ)も使います。
kramdownは、markdownの原稿をhtml, LaTeXに変換する機能を提供するものです。これを利用して、同じ内容のレポートをhtml, LaTeXの両形式で作成することを考えます。
LaTeX文書は、dvipdfmxコマンドを用いて、最終的にはpdfにすることを念頭に置いて作成します。
圧縮ファイル lm_s.zip に含まれているサンプルスクリプト lm_htm01.rb などを私が動かした環境は次のとおり。
素材データ data01.csv がカレントディレクトリにある状態で、次のように実行します。
ruby lm_htm01.rb [enter]
そうすると、result201020.htm, cor201020.png, lm201020.png の３つのファイルが書き出されます。後者の２つはグラフの画像ファイルです。
result201020.htm をブラウザで開くと、レポートを閲覧できます。
「201020」の数字列は、スクリプトを実行した時の時・分・秒を示します。その時々でかわります。
ある喫茶店でのアイスティーの売上げと最高気温が列記されたcsvデータがあるとします。１日ごとにデータが記録されています。
最高気温が高いほどアイスティーの売上げも上がる、そんな関係を予想したとして、この２つをどのような関係式で示せるかを検証します。
実は、この事例は、下のサイトで詳しく解説されています。
取り上げられているデータは、「マンガでわかる統計学 回帰分析編」のデータのようです。
圧縮ファイルに同梱してある data01.csv は、基本的にそれと同じです。
上記サイトを私なりに咀嚼したつもりで(？)、簡単なレポートにしてみました。そのレポートは次のwebで閲覧できます。同文のpdfもあります。
pdf版：最高気温とアイスティーの売上げの関係に関する回帰分析
上のwebを作成するためのrubyスクリプトについて、以下で記します。
説明の主眼は、Rの分析結果をどのようにして取り込むか、また、それをhtml, LaTeXの文書にどのようにして変換するかというところにあります。回帰分析の意味や方法については、前述のサイトを参照して下さい。
Rプログラムでは、データの読み込み、相関のチェック、回帰分析の結果取得、それと散布図等のグラフ作成を行います。グラフはpngファイルを書き出します。
相関のチェックは cor.test()、回帰分析は lm() で行います。
それら材料を受け取って、実際のレポートにするのは ruby の役割です。
下にRプログラムを掲げます。
−−−− Rプログラムここから
dtf <- read.csv("ddd", header=T)
res <- cor.test(dtf$xxx, dtf$yyy) # 相関の検定
robj(res, "相関")
# 散布図の作成
png(file="cor<%= id_n %>.png")
par(family = "Japan1GothicBBB") # 日本語フォントにゴチックを指定
par(mar = c(5, 5, 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="lm<%= id_n %>.png")
par(family = "Japan1GothicBBB") # 日本語フォントにゴチックを指定
par(mar = c(5, 5, 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()
−−−− Rプログラムここまで
上のRプログラムの中に ddd, xxx, yyy のように、小文字のアルファベットが３つ重なるものが出てきます。
これは、後でrubyの正規表現置換によって「<%= ddd %>」とか「<%= xxx %>」などに置換されます。つまりerbで処理できる形に改められるわけです。
最初からerb対応の形で記述しておくのが正当なやり方だと思いますが、その簡便さに捨てがたいものがあり、こんなふうにしました。場合によってはトラブルの原因になるかもしれないので、同じやり方を採る時は注意して下さい。
予想がつくと思いますが、最終的には、dddが読み込むデータファイル名(data01.csv)、xxxが説明変数(最高気温)、yyyは目的変数(アイスティーの売上げ)に置き換えられます。
「<%= id_n %>」というのが出てきますが、id_nは識別番号です。
グラフの画像ファイル名を単に cor.png とか lm.png としてもいいのですが、cor01.png, lm01.png のようにした方が、別のレポート作成と混同せずに済むと思って付けたものです。
実際のサンプルスクリプトでは、この id_n に時刻の「時・分・秒」を割り当てます。スクリプトを実行した時刻です。８時19分32秒であれば 081932 となります。
グラフは２つ書き出します。図１が散布図、図２は散布図に回帰直線を書き入れたものです。
相関のチェックには cor.test() を使いましたが、回帰分析のところで決定係数とその有意性を見るので、cor() で相関係数の値をチェックするので十分かもしれません。
というより、説明変数が２つ以上になる重回帰分析の場合は、それらすべての組合せについて相関を一覧することのできる cor() を用いるのがいいと思います(重回帰分析の例参照)。
回帰分析の結果については、lm() の戻り値と、それをsummary() にかけたものの２つをrubyに引きわたしています。
summary() にかけたものだけでもいいのですが、lm() の戻り値の方が、回帰係数を容易に取り出せるので、それも引きわたすようにしました。
Rのオブジェクトをrubyに引きわたすのは、robj() で行っています。
なお、lm() の戻り値や、それを summary() にかけたものには、決定係数の有意性を見るためのp値が盛り込まれていません。
そこで、そのp値を算出して変数resに付け加えます。その上でresをrubyに引きわたしています。
rubyスクリプトの大まかな流れは次のようになります。
上の第１ステップ：data_name, xxx, yyy の設定を変更すれば、他の部分を変更しなくても、別のデータに関する回帰分析のレポートを作成できます。xxxとyyyが高い相関関係にあるとの前提で書かれたレポートなので、内容的に的外れになってしまう可能性はありますが、ともあれ同じ雛形のレポートが作成されます。
圧縮ファイルに同梱してある lm_htm02.rb は、web上のExcelファイルを材料にして、BMIと体脂肪率の関係を分析するものです。
また、最後の「kramdownによるhtmlの生成」のところを少し書き換えれば、htmlでなくLaTeX文書の作成になります。
同梱の lm_tex01.rb を実行すると、同じ内容のLaTeX文書が書き出されます。
必要なライブラリをrequireした後に次の３行があります。
data_name = "data01.csv" # 素材データ
xxx = "最高気温" # 説明変数(推定の材料)
yyy = "アイスティーの売上げ" # 目的変数(推定したい項目)
これらが回帰分析を行う上での基本情報です。
この３行を下のように書き換えると lm_htm02.rb になります。
data_name = "http://www.qmss.jp/databank/preliminaries/data/1-2b-2.xls"
xxx = "BMI"
yyy = "体脂肪率"
上記のExcelファイルは、文科省基準「新体力テスト」（12−19歳）のデータのようです。
「わかりやすい統計学」で取り上げているもので、多変量データの解説で参考にしているデータです。重回帰分析にも使えると思います。
この図書の紹介サイトは次のところです。様々なExcelデータをここからダウンロードできます。
Rプログラムで読み込むファイルの名前を変数 ddd に代入します。
csvファイルであれば素直に「ddd = "data01.csv"」のようにできますが、Excelファイルの場合は、素材データが書かれているワークシートをcsvファイルとして書き出して(テンポラリファイルに書き出す)、そのファイル名をdddに代入します。
該当部分は次のとおり。
if data_name =~ /\.xls$/i
sheets = Rrx.xls2ary(data_name) # Excelデータを配列に変換
tcsv = Rrx.temp_make(sheets[0], ".csv")
ddd = tcsv.path
else
ddd = data_name
end
「tcsv = Rrx.temp_make(sheets[0], ".csv")」は、第１ワークシートをcsvファイルとして書き出すための記述です。
そのcsvファイルのフルパス名は tcsv.path で参照できます。
tcsv.unlink とすれば、csvファイルを削除できます。Rプログラムを実行した後であれば、削除してかまいません。
上の記述では拡張子 *.xls のExcelファイルだけを対象にしていますが、rubyライブラリのrooをrequireしていれば、*.xlsx, *.ods なども扱えます。
MS-Windows環境下でExcelがインストールされているなら、拙作exlap(ruby用のライブラリ)をrequireすると、同じように *.xlsx, *.xlsm, *.ods などを扱うことができます。
*.xls を扱うだけなら spreadsheet をrequireします。MS-Windows以外の環境でも使えます。
Rプログラムとmarkdown原稿は、rubyスクリプトの最後の方に書いてあります。具体的には「__END__」という行の後ろに書いてあります。なので、DATA.read で読み込むことができます。
Rプログラムとmarkdown原稿の区切りは、「###### 以下はmarkdownの原稿」という行です。行頭に半角のシャープ記号が３個以上あると、その行を区切りとみなして分割します。
分割する前に、ddd, xxx, yyy などの小文字アルファベット３段重ねを「<%= ddd %>」などに置換します。この置換は、Rプログラムとmarkdown原稿に共通するので、分割前に行います。
アルファベット３文字に続けて半角数字がくる場合(xxx2など)も置換の対象になります。また、３文字の後に半角下線と英数字がくるケース(xxx_mainなど)、３文字の後に「[……]」がくるケース(xxx["分析結果"]など)も置換の対象になります。
該当の処理部分は下のとおり。
src_data = DATA.read # Rプログラムとmarkdownの雛形文書
reg = /([a-z])\1\1(\d+|_[A-Za-z0-9]+|[\[].+?[\]])*/
src_data = src_data.gsub(reg, "<%= \\& %>") # aaa, xxx などの変換
r_prog, mkd_str = src_data.split(/^###+.*?\n+/) # Rとmarkdownの分離
id_n = Time.now.strftime("%H%M%S") # 時・分・秒:出力ァイル名に盛り込む
Rプログラムが変数r_prog、markdown原稿が変数mkd_strにセットされます。
markdown原稿の先頭にある空白行は消去されます。
r_progをerbにかけて最終的に整え、その上で Rrx.rexec(r_prog) で実行します。
実行結果は変数 hs にセットします。hsはハッシュ(連想配列)です。
一応、該当箇所を掲げておきます。
r_prog = ERB.new(r_prog).result(binding)
hs = Rrx.rexec(r_prog) # Rの実行
tcsv.unlink if tcsv # csvテンポラリファイルの削除
Rプログラムが実行中に警告メッセージなどを出すかもしれない、と思った場合は、次の１行を挿入します。
print hs[:log] こうすると、Rプログラムの実行結果がすべて標準出力に出力されます。警告メッセージなども出力されます。
変数hsを参照して、Rの分析結果をrubyの適当な変数にセットします。
相関の検定結果は cor に代入します、その中でも、相関係数・自由度・p値を cor_result にセットします。
回帰分析の結果は lm、回帰分析の要約(summary)は lm2 にセットします。
この lm, lm2 を材料にして、y切片をaaa、xの傾きをbbb、決定係数をrrr、回帰係数に関する一連の情報(２次元配列)をcoefにセットします。
coefは、２次元配列なので表形式で出力するのが容易ですが、項目名が英語です。例えば、標準誤差が「Std. Error」となっています。これを日本語に変更します。
該当箇所は次のとおり。
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] = "定数項"
「定数項」は intercept を訳したものですが、y切片のことです。
まず、変数 mkd_str (markdownの原稿)をerbにかけて、最終調整します。
次に、htmlのスタンドアロンを生成するためのテンプレートをテンポラリファイルとして書き出します。「Rrx.html_template()」で行えます。
あとは kramdown で変換処理を行い、結果をファイルに出力します。
該当部分は次のとおり。
mkd_str = ERB.new(mkd_str).result(binding)
tfile = Rrx.html_template() # テンプレート生成(tempfile)
kdoc = Kramdown::Document.new(mkd_str, :template=>tfile.path)
html_str = kdoc.to_html
File.open("result#{id_n}.htm", "w") {|ff| ff.write html_str}
tfile.unlink # テンプレートファイルを削除
この部分を適当に変更すると、LaTeX文書を出力するものになります。
lm_tex01.rb の該当箇所を掲げてみます。
mkd_str = ERB.new(mkd_str).result(binding)
tfile = Rrx.tex_template() # テンプレート生成(tempfile)
kdoc = Kramdown::Document.new(mkd_str, :template=>tfile.path)
tex_str = kdoc.to_latex
tex_str = tex_str.sub(/^(\\usepackage)(\{graphicx\})/, "\\1[dvipdfmx]\\2")
File.open("result#{id_n}.tex", "w") {|ff| ff.write tex_str}
tfile.unlink # テンプレートファイルを削除
テンプレート生成のところが html_template() から tex_template() になっているのと、変換指定の to_html が to_latex になっている点が相違点です。
その他、変換結果の tex_str に対して sub() で置換処理を施していますが、これは dvipdfmx コマンドでpdfを生成する時に、画像をうまく扱えるようにするためのものです。
\usepackage{graphicx}
上の１行を少しだけ書き換えて、次のようにします。
\usepackage[dvipdfmx]{graphicx}
こうすると、最終的に生成されるpdfに画像がうまく取り込まれます。
ちなみに、lm_tex01.rb では、グラフの画像ファイルをpngではなくpdfで生成するようにしています。そのため、Rプログラムとmarkdownの原稿部分のpngの３文字がpdfに変更されています。
pngのままでもいいかもしれませんが、pdfの方がスムーズに処理されると思います。
LaTeX文書を基にしてpdfファイルを生成する方法については、次のサイトで紹介しているので、必要に応じて参照して下さい。
kramdownに関する覚え書き 〜 markdown, html, LaTeX の変換処理
以上がrubyによる処理です。
html用の原稿とLaTeX用の原稿には僅かな違いがあります。先にもちょっと触れましたが、html用ではグラフをpng形式で書き出し、LaTeXではpdf形式で書き出します。
しかし、違いはそれだけです。html用の中のpngの３文字をpdfに置換すると、LaTeX用と同じものになります。kramdownによって変換処理するため、ほぼ同じ原稿から html, LaTeX の２種類を生成できます。
もちろん、実際には個別の調整を加えたくなる点が出てきます。例えば、今回生成されるpdfでは、グラフの画像が文章中に挿入されるのではなく、まるまる１ページを占領します。これを避けようと思うと、画像の大きさを調整したり、「なるべく文中に挿入するよう努めなさい」という命令をLaTeX文書に埋め込む必要があります。
とはいえ、基本的に同じ原稿から html, LaTeX の２種類を作ることができるのは、いろいろと便利です。
markdownの原稿をこの後に掲げますが、その中に、普通のmarkdown記述にはないものとして次の１行があります。
<%= Rrx.matrix2table(coef, 'border="1"') %>
matrix2table() は、Rから引きわたされた２次元配列(ここでは coef)をmarkdownの表に変換するメソッドです。markdownの表は、例えば次のようなものです。
|名称|電話番号|
|時刻|117|
|天気予報|177|
サンプルに出てくるcoefには、回帰係数の値、標準誤差、t検定量、t検定のp値(両側)が２次元配列としてセットされています。これを表形式で出力します。
今回は単回帰分析なので、見出しを除くと２行だけの表ですが、説明変数が２つ以上になると(つまり重回帰分析の場合)、この表の行数が増えます。
markdown記法については説明を省略します。多くの解説サイトがあるので、必要に応じて参照して下さい。
以下、長くなりますがmarkdownの原稿を掲げておきます。このサイトは、それでおわりです。
−−−−−−−− 原稿ここから
# xxxとyyyの関係に関する回帰分析
ここにxxxとyyyを列記したデータがある。これについて回帰分析を行う。
xxxを材料にしてyyyを推定できるか調べる。
このような推定を専門的表現で「yyyをxxxに回帰させる」というらしい。
## 相関の確認
まず、回帰分析を行う意味のあるデータなのかどうかを確認する。
２つの項目の相関が認められなければ、
回帰分析の結果が実質的な意味を持たない。
「xxx」と「yyy」の相関係数を確認すると次のとおり。
<%= cor_result %>
有意な強い相関が認められる。
ちなみに、散布図は図１のようになる。
* * * * * * * *
## 回帰分析
強い相関が認められるので、回帰分析を行ってみた。その結果は次のとおり。
- y切片： aaa
- xの傾き： bbb
つまり、次のような推定のための計算式(回帰式)が求められる。
「yyy」＝aaa＋bbb×「xxx」
上の回帰式によって描かれる直線が回帰直線である。
散布図に回帰直線を書き入れてみると図２のようになる。
個々の実測値と回帰式から得られる値との間には誤差(残差)がある。
すべての実測値についてこの誤差の平方和が最小になるように
考え出されたのが上の回帰式である。
* * * * * * * *
## 回帰式の信頼性：決定係数
先の回帰式がどの程度信頼できるかをみる一つの手がかりは、決定係数である。
- 決定係数: rrr
- Ｆ検定の値: <%= lm2["fstatistic"][0] %> (決定係数に関するＦ検定量)
- Ｆ検定のp値: <%= lm2["f_p.value"] %>
決定係数は、１に近いほど信頼性が高い。
今回の rrr という値は、「yyy」の変動が「xxx」の変動によって
rrr_main％だけ説明できることを意味している。
残りのrrr_other％は、回帰式に含まれていない別の要素が関係しているとみられる。
決定係数の信頼性はＦ検定の結果により判断する。もっともっと幅広く
調査したとき、もしかすると決定係数が０になるかもしれない。その危険性を
示すのがＦ検定のp値である。これが0.05(5%)とか0.01、あるいは0.001より
小さければ、そのレベルに応じた有意性が認められることになる。
このp値が0.05以上だと、「母集団において決定係数が０である」
という仮説を、母集団に関して棄却できないことになる。
## 回帰式の信頼性：回帰係数
回帰係数(coefficient：y切片やxの傾きの関連)
に関する詳細な情報を掲げると次のとおり。
<%= Rrx.matrix2table(coef, 'border="1"') %>
この表で、１行目の「定数項」の「見積り」はy切片の値であり、
２行目の「<%= coef[2][0] %>」の「見積り」の欄はxの傾きの値である。
重回帰分析の場合は、説明変数の個数に応じて、
３行目・４行目・……が表示されることになる。
表の右端の「t検定時のp値」は、該当の変数が実は影響力を
持たないかもしれない危険性(確率)を示す。
xの傾きが０であれば、xの値が何であってもyの値に影響しない。
xxxの傾きがbbbと見積もられてはいるが、もっともっと調査対象を拡げると、
実は０ということもあるかもしれない。その「かもしれない」の確率がp値である。
p値が十分小さければ、有意性が認められることになる。
定数項のp値については、実際に言及することは少ないと思うが、
もっともっと幅広く調査した時にy切片が０であるかもしれない確率を示す。
−−−−−−−− 原稿ここまで
〜 以上 〜
Copyright (C) T. Yoshiizumi, 2013 All rights reserved.