Recoding a factor variable in a data.frame by the levels() function
When you want to change the coding of a factor variable in R, you can use the levels()
function.
Recoding z
as 1 = A, 2-6 = B, 7 = C
d <- data.frame(z = gl(7, 1, 14)) d$zc <- factor(d$z) levels(d$zc) <- list(A = 1, B = 2:6, C = 7) d
Result
> levels(d$zc) <- list(A = 1, B = 2:6, C = 7) > d z zc 1 1 A 2 2 B 3 3 B 4 4 B 5 5 B 6 6 B 7 7 C 8 1 A 9 2 B 10 3 B 11 4 B 12 5 B 13 6 B 14 7 C
Recoding z
as 1 = "1", 2-6 = "2-6", 7 = "7"
levels(d$zc) <- list("1" = 1, "2-6" = 2:6, "7" = 7) d
Result
> levels(d$zc) <- list("1" = 1, "2-6" = 2:6, "7" = 7) > d z zc 1 1 1 2 2 2-6 3 3 2-6 4 4 2-6 5 5 2-6 6 6 2-6 7 7 7 8 1 1 9 2 2-6 10 3 2-6 11 4 2-6 12 5 2-6 13 6 2-6 14 7 7
Penalized quasi-likelihood について
GLMM の推定アルゴリズムに Penalized quasi-likelihood (PQL) という方法があるのですが、整理のためメモを書くことにしました。
PQL
混合効果モデルのパラメータ推定において、周辺尤度の積分を必要とする方法があることはよく知られていると思います。
GLMM ではリンク関数が非線形なものを扱うが故に、この積分が非常に大変でした (現代の計算機では苦労せず解けるものも多い)。
一昔前はこの積分が大変だったため、非線形な関数を線形な関数に近似してあげて、LMM の範疇で扱えるようにした方法が PQL です。
GLMM として、確率変数 $\mathbf{Y}_i$ の期待値が、変量効果 $\mathbf{b}_i$ を与えた下で、
\[
\mathrm{E}[\mathbf{Y}_i|\mathbf{b}_i]=g^{-1}(\mathbf{X}_i\boldsymbol{\beta}+\mathbf{Z}_i\textbf{b}_i)=g^{-1}(\boldsymbol{\xi}_i)=\boldsymbol{\mu}_i
\] とモデル化された一般的なものを考えます ($\mathbf{X}_i$, $\mathbf{Z}_i$ はデザイン行列、$\boldsymbol{\beta}$ は固定効果のパラメータベクトル、$\boldsymbol{\xi}_i = \mathbf{X}_i\boldsymbol{\beta}+\mathbf{Z}_i\textbf{b}_i$、$g^{-1}$ はリンク関数の逆関数)。
PQL では $\boldsymbol{\xi}_i$ の関数である $\boldsymbol{\mu}_i$ に対して $\boldsymbol{\beta}=\hat{\boldsymbol{\beta}}$, $\mathbf{b}_i=\hat{\mathbf{b}_i}$ の周りの一次の Taylor 展開、
\[
g^{-1}(\mathbf{X}_i\boldsymbol{\beta}+\mathbf{Z}_i\textbf{b}_i) \approx g^{-1}(\mathbf{X}_i\hat{\boldsymbol{\beta}}+\mathbf{Z}_i\hat{\textbf{b}_i}) + \hat{\Delta}_i\mathbf{X}_i(\boldsymbol{\beta}-\hat{\boldsymbol{\beta}}) + \hat{\Delta}_i\mathbf{Z}_i(\textbf{b}_i-\hat{\textbf{b}_i})
\] を考えます、ただし、
\[
\hat{\Delta}_i= \mathrm{diag} \left( \frac{\partial g^{-1}(\xi_{ij})}{\partial \xi_{ij}} \Big|_{\boldsymbol{\xi}_i=\hat{\boldsymbol{\xi}}_i}\right)_{1 \leq j \leq n_i}
\] とします ($n_i \times n_i$ 次元行列)。
これを整理すると、
\[
\hat{\Delta}_i^{-1}(\boldsymbol{\mu}-g^{-1}(\mathbf{X}_i\hat{\boldsymbol{\beta}}+\mathbf{Z}_i\hat{\textbf{b}_i}))+ \mathbf{X}_i\hat{\boldsymbol{\beta}}+ \mathbf{Z}_i\hat{\textbf{b}_i} \approx \mathbf{X}_i\boldsymbol{\beta}+ \mathbf{Z}_i\textbf{b}_i
\] が得られます。
左辺をよく見ると、これは、$\hat{\Delta}_i^{-1}(\mathbf{Y}_i-g^{-1}(\mathbf{X}_i\hat{\boldsymbol{\beta}}+\mathbf{Z}_i\hat{\textbf{b}_i}))+ \mathbf{X}_i\hat{\boldsymbol{\beta}}+ \mathbf{Z}_i\hat{\textbf{b}_i} \equiv \mathbf{P}_i$ の期待値になっていることが分かります。
確率変数 $\mathbf{Y}_i$ 以外は全て定数なので、新しい作業確率変数として $\mathbf{P}_i$ を定義します。
また、同じく $\mathbf{Y}_i$ 以外が全て定数なので、$\mathrm{Var}(\mathbf{P}_i\mid\mathbf{b}_i)= \hat{\Delta}_i^{-1} \mathrm{Var}(\mathbf{Y}_i\mid\mathbf{b}_i) \hat{\Delta}_i^{-1}$ が簡単に分かります。
そして、右辺からは非線形なリンク関数がうまく消去できていて、線形化することができました。
以上の近似を使って、LMM、
\[
\mathbf{P}_i=\mathbf{X}_i\boldsymbol{\beta}+\mathbf{Z}_i\mathbf{b}_i+\boldsymbol{\epsilon}^*_i,~
\mathrm{Var}(\boldsymbol{\epsilon}^*_i)= \mathrm{Var}(\mathbf{P}_i \mid \textbf{b}_i)
\] を当てはめたものが、PQL と呼ばれる方法です。
また、類似の Marginal quasi-likelihood (MQL) という方法では、変量効果の方を母集団平均で置き換えるため、$\mathbf{b}_i=\mathbf{0}$ を用います。
ただし、MQL は変量効果が非常に小さい場合でないとうまくいかないことが知られています。
また、アウトカムが二項分布に従う場合は、PQL と MQL による近似は特にうまくいかないケースがあることが知られています。
感想
とまあ、一次の Taylor 展開を使ったシンプルな方法だという事が分かります。
計算は早くなりますが、デメリットもあります。
作業確率変数 $\mathbf{P}_i$ に置き換えるため、$\mathbf{y}_i$ に基づく尤度が定義できない事はデメリットの一つでしょう。
どこかに quasi-likelihood を使う方法だから良くないと書いてあった気がしますが、そうではないんです。
本質的な問題点は作業確率変数を使った近似の方ですよね。
よって、利用可能であれば、適応型ガウス–エルミート求積法などの周辺尤度を直接計算するような方法が望ましいでしょう。
ちなみに SAS 9.3 の GLIMMIX Procedure は METHOD = RSPL (PQL で REML 推定) がデフォルトです。
参考文献
[1] Breslow N. Whither PQL? UW Biostatistics Working Paper Series. Working Paper 192. http://biostats.bepress.com/uwbiostat/paper192
Log-gamma distribution
A log-gamma random variable, $Y$, is defined as,
\[
Y=\exp(-X)
\] where $X \in [0, \infty)$ follows a gamma distribution, and $Y \in (0, 1]$.
The probability density function of the log-gamma distribution is
\[
f_Y(y)=\frac{b^a\left\{\log(1/y) \right\}^{a-1}y^{b-1}}{\Gamma(a)},
\] where $b > 0$, and $a > 0$.
The expected value of $Y$ is
\[
\mathrm{E}[Y]=\int_0^1 yf_Y(y) \mathrm{d}y=\left( \frac{b}{1+b} \right )^a.
\] The variance of $Y$ is
\[
\mathrm{Var}[Y]=\mathrm{E}[Y^2]-(\mathrm{E}[Y])^2= \left( \frac{b}{2+b} \right)^a - \left( \frac{b}{1+b} \right)^{2a}.
\]
R code:
# shape = a, rate = b a <- 4 b <- 7 r <- -rgamma(100000, shape = a, rate = b) mean(exp(r)) (b/(1 + b))^a mean(exp(r - log((b/(1 + b))^a))) var(exp(r)) (b/(2 + b))^a - (b/(1 + b))^(2*a)
Result:
> mean(exp(r)) [1] 0.5855565 > (b/(1 + b))^a [1] 0.5861816 > > mean(exp(r - log((b/(1 + b))^a))) [1] 0.9989336 > > var(exp(r)) [1] 0.02245031 > (b/(2 + b))^a - (b/(1 + b))^(2*a) [1] 0.0223414
Note
\[
\int_0^1 yf_Y(y) \mathrm{d}y= \frac{b^a}{\Gamma(a)} \int_0^1 \left\{\log(1/y) \right\}^{a-1}y^b \mathrm{d}y
\] \[
\begin{align*} z &= \log(1/y) \\
\mathrm{d}z &= -\frac{1}{y} \mathrm{d}y \\
\Re(z) &> 0, \log(1/0)=\infty, \log(1/1)=0.
\end{align*}
\] \[
\begin{align*}
\int_0^1 yf_Y(y) \mathrm{d}y &= \frac{b^a}{\Gamma(a)} \int_{\infty}^0 -z^{a-1} e^{-z(1+b)} \mathrm{d}z \\
&= \frac{b^a}{\Gamma(a)} \frac{\Gamma(a)}{(1+b)^a} \\
&= \left( \frac{b}{1+b} \right)^a \end{align*}.
\]
Java Compiler API を使って何かしてみよう
Java のバイナリから Java ソースコードをコンパイルするための API を使って、String 型の文字列をコンパイルし、クラスから static なメソッドを呼んで値を返す方法を模索していました。
以下のプログラムと解説を参考にしています
- Compiling from Memory : Java Compiler tools « JDK 6 « Java
- Javaリフレクションメモ(Hishidama's Java Reflection Memo)
import java.io.IOException; import java.lang.reflect.InvocationTargetException; import java.lang.reflect.Method; import java.net.URI; import java.util.Arrays; import javax.tools.Diagnostic; import javax.tools.DiagnosticCollector; import javax.tools.JavaCompiler; import javax.tools.JavaFileObject; import javax.tools.SimpleJavaFileObject; import javax.tools.ToolProvider; import javax.tools.JavaCompiler.CompilationTask; public class CompileSourceInMemory { public static void main(String args[]) throws IOException { JavaCompiler compiler = ToolProvider.getSystemJavaCompiler(); DiagnosticCollector<JavaFileObject> diagnostics = new DiagnosticCollector<JavaFileObject>(); String src = "public class HelloWorld {" + "public static int main(String args[]) {" + "System.out.println(\"This is in another java file\");" + "return(1);" + "}" + "}"; JavaFileObject file = new JavaSourceFromString("HelloWorld", src); String[] compileOptions = new String[] {"-d", "bin"} ; Iterable<String> compilationOptionss = Arrays.asList(compileOptions); Iterable <? extends JavaFileObject > compilationUnits = Arrays.asList(file); CompilationTask task = compiler.getTask( null, null, diagnostics, compilationOptionss, null, compilationUnits ); boolean success = task.call(); for (Diagnostic<?> diagnostic : diagnostics.getDiagnostics()) { System.out.println(diagnostic.getCode()); System.out.println(diagnostic.getKind()); System.out.println(diagnostic.getPosition()); System.out.println(diagnostic.getStartPosition()); System.out.println(diagnostic.getEndPosition()); System.out.println(diagnostic.getSource()); System.out.println(diagnostic.getMessage(null)); } System.out.println("Success: " + success); if (success) { Object ret; try { Class<?> clazz = Class.forName("HelloWorld"); Method method = clazz.getMethod("main", new Class[] { String[].class }); ret = method.invoke(null, new Object[] { null }); System.out.println(ret); } catch (ClassNotFoundException e) { System.err.println("Class not found: " + e); } catch (NoSuchMethodException e) { System.err.println("No such method: " + e); } catch (IllegalAccessException e) { System.err.println("Illegal access: " + e); } catch (InvocationTargetException e) { System.err.println("Invocation target: " + e); } } } } class JavaSourceFromString extends SimpleJavaFileObject { final String code; JavaSourceFromString(String name, String code) { super(URI.create("string:///" + name.replace('.', '/') + Kind.SOURCE.extension), Kind.SOURCE); this.code = code; } @Override public CharSequence getCharContent(boolean ignoreEncodingErrors) { return code; } }
ポイント
私の環境では Eclipse 上で実行しているので、
String[] compileOptions = new String[] {"-d", "bin"} ;
を入れています。
もう一つは、
Object ret; try { Class<?> clazz = Class.forName("HelloWorld"); Method method = clazz.getMethod("main", new Class[] { String[].class }); ret = method.invoke(null, new Object[] { null }); System.out.println(ret); }
という感じにリフレクションを使って、返り値がある場合に拡張してみました。
とりあえず Eclipse 上では実行できた。ので、jar ファイルなどにしても動けばいいのですが・・・
2012/06/25 追記
実行可能 jar ファイル化して動かなかった原因は、jdk ではないと動かない (tools.jar 必須) のに、jre から動いていたため・・・
最新 jdk 以外を全てアンインストール / jar ファイルの関連づけを変更したら普通に動きました。
とりあえず、
JavaCompiler compiler = ToolProvider.getSystemJavaCompiler(); if (compiler == null) { System.out.println("null dayo!"); }
で、compiler を取ってこれているかどうかはチェックはしておいた方が良いでしょう。
null の場合は tools.jar が読み込めていないので、jdk から実行できていそうかなどをちゃんとチェックしておくとよいみたいです。
2012/06/26 リフレクションメモ
複数引数のリフレクションサンプルを作ってみました。
可変長引数がないので、以下の様に列挙するだけで簡単に使えますね。
import java.lang.reflect.InvocationTargetException; import java.lang.reflect.Method; public class a { public static void main(String[] args) { Class<?> clazz; try { clazz = Class.forName("test"); Object obj = clazz.newInstance(); Method method = obj.getClass().getMethod("fx", double.class, double.class, double[].class); System.out.println((Double) method.invoke(null, 2.3, 1.0, new double[] {1.0, 2.0, 1.0})); // 7.3 が出力される } catch (ClassNotFoundException e) { e.printStackTrace(); } catch (InstantiationException e) { e.printStackTrace(); } catch (IllegalAccessException e) { e.printStackTrace(); } catch (SecurityException e) { e.printStackTrace(); } catch (NoSuchMethodException e) { e.printStackTrace(); } catch (IllegalArgumentException e) { e.printStackTrace(); } catch (InvocationTargetException e) { e.printStackTrace(); } } } class test { public static double fx(double a, double b, double[] c) { double cc = 0; for (int i = 0; i < c.length; i++) { cc += c[i]; } return a + b + cc; } }
R Commander による Kaplan-Meier プロットの作図手順
- R のインストール バージョン 2.15.0 以上の R をインストールします。ソフトウェアのインストールをした経験がある方なら特に困ることはないでしょう。
- Windows の場合
- Mac の場合
- http://www.okada.jp.org/RWiki/?R%20%A4%CE%A5%A4%A5%F3%A5%B9%A5%C8%A1%BC%A5%EB に詳しい情報がありますので、困ったら参照しましょう
- R Commander のインストール
- R Commander プラグインのインストール
- R Commander の起動
- R Commander プラグインの読み込み
- プラグインのロード
- 再起動を促されるので再起動します
- データの読み込み
- イベント変数をイベントありの場合が 1、打ち切りの場合が 0 でコーディングした形式のデータを準備し、クリップボードにコピーしておきます
- データ、データのインポート、テキストファイル・・・を選択します
- データセット名を指定し、変数名はあり、クリップボード、タブをチェックし、OK をクリックします
- 作図
R を起動し、以下のコマンドを実行するか、パッケージのインストール画面等から "Rcmdr" を選択します。
install.packages("Rcmdr")
インストールエラーが表示された場合は、もう一度同じ操作を行ってみてください。
同様にプラグインをインストールします。
以下のコマンドを実行するか、パッケージのインストール画面等から "RcmdrPlugin.KMggplot2" を選択します。上の画像のように、複数のパッケージを同時にインストールすることもできます。
install.packages("RcmdrPlugin.KMggplot2")
以下のコマンドを実行するか、パッケージのロード画面から "Rcmdr" を選択します。
library("Rcmdr")
RcmdrPlugin.KMggplot2 0.2-0 以降では、R Commander の起動とプラグインのロードを同時に行うことができます。
以下のコマンドを実行するか、パッケージのロード画面から "RcmdrPlugin.KMggplot2" を選択します。
library("RcmdrPlugin.KMggplot2")
ロードが終了したら、5 番目のステップ (プラグインの読み込み) を実行する必要はありません。
グラフ作成に必要なプラグインを読み込みます。







記号微分
deriv
を利用して、R で記号微分してみます。
正規分布の対数尤度関数を記号微分してみよう
正規分布の対数尤度関数は、
\[
l(\mu, \s; x)=-\frac{\log(2\pi\s^2)}{2}-\frac{(x - \mu)^2}{2\s^2}
\] ですので、これを偏微分すると、
\[
\begin{align*}
\frac{\partial l(\mu, \s; x)}{\partial \mu} &= \frac{x-\mu}{\s^2} \\
\frac{\partial l(\mu, \s; x)}{\partial \s} &= -\frac{1}{\s}+\frac{(x-\mu)^2}{\s^3}
\end{align*}
\] となる予定です。
deriv
を用いた数値微分では、まず被微分関数を expression
、call
、または formula
の形で定義します。
fx <- expression( -log(2*pi*sd^2)/2 - (x - mean)^2/(2*sd^2) )
あとは、
を指定して実行するだけです。
dfx <- deriv( fx, c("mean", "sd"), function.arg = c("x", "mean", "sd") ) dfx
結果は以下のようになりました。
> dfx function (x, mean, sd) { .expr1 <- 2 * pi .expr2 <- sd^2 .expr3 <- .expr1 * .expr2 .expr7 <- x - mean .expr8 <- .expr7^2 .expr9 <- 2 * .expr2 .expr14 <- 2 * sd .value <- -log(.expr3)/2 - .expr8/.expr9 .grad <- array(0, c(length(.value), 2L), list(NULL, c("mean", "sd"))) .grad[, "mean"] <- 2 * .expr7/.expr9 .grad[, "sd"] <- -(.expr1 * .expr14/.expr3/2 - .expr8 * (2 * .expr14)/.expr9^2) attr(.value, "gradient") <- .grad .value }
一応確認しておくと、
.expr2 <- sd^2 .expr7 <- x - mean .expr9 <- 2 * .expr2 .grad[, "mean"] <- 2 * .expr7/.expr9 # 2 * .expr7/.expr9 = # 2 * (x - mean) / 2 * sd^2 = # (x - mean) / sd^2
および、
.expr1 <- 2 * pi .expr2 <- sd^2 .expr3 <- .expr1 * .expr2 .expr7 <- x - mean .expr8 <- .expr7^2 .expr9 <- 2 * .expr2 .expr14 <- 2 * sd .grad[, "sd"] <- -(.expr1 * .expr14/.expr3/2 - .expr8 * (2 * .expr14)/.expr9^2) # -.expr1 * .expr14/.expr3/2 = # -2 * pi * 2 * sd / (2 * pi * sd^2) / 2 = # -1/sd # .expr8 * (2 * .expr14)/.expr9^2 = # (x - mean)^2 * (2 * 2 * sd) / (2 * sd^2)^2 = # (x - mean)^2 / sd^3
となります。
このぐらいなら手計算の方が早いですが、いろいろな場面で利用できるのではないでしょうか。
また、Hessian を自動で出してくれる deriv3
もあります。
注意
deriv
のヘルプのタイトルは、
Symbolic and Algorithmic Derivatives of Simple Expressions
となっていて、以下の説明があります。
The internal code knows about the arithmetic operators +, -, *, / and ^, and the single-variable functions exp, log, sin, cos, tan, sinh, cosh, sqrt, pnorm, dnorm, asin, acos, atan, gamma, lgamma, digamma and trigamma, as well as psigamma for one or two arguments (but derivative only with respect to the first). (Note that only the standard normal distribution is considered.)
RcmdrPlugin.KMggplot2_0.1-0 is on CRAN now
I posted a new version of the "RcmdrPlugin.KMggplot2" package, which is an Rcmdr plug-in for a "ggplot2" GUI front-end.
This package assists you to make "ggplot2" graphics.
- RcmdrPlugin.KMggplot2 (CRAN)
NEWS
Changes in version 0.1-0 (2012-05-18)
- Restructuring implementation approach for efficient maintenance.
- Added
options()
for storing package specific options (e.g., font size, font family, ...). - Added a theme:
theme_simple()
. - Added a theme element:
theme_rect2()
. - Added a list box for
facet_xx()
functions in some menus (Thanks to Professor Murtaza Haider). - Kaplan-Meier plot: added confidence intervals.
- Box plot: added violin plots.
- Bar chart for discrete variables: deleted dynamite plots.
- Bar chart for discrete variables: added stacked bar charts.
- Scatter plot matrix: added univariate plots at diagonal positions (ggplot2::plotmatrix).
- Deleted the dummy data for histograms, which is large in size.