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

2012-09-01

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

2012-07-26

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*}.
\]

Reference

[1] Saberali SA, Beaulieu NC. Calculating the distribution of sums of log-gamma random variables.

2012-06-23

Java Compiler API を使って何かしてみよう

Java

Java のバイナリから Java ソースコードをコンパイルするための API を使って、String 型の文字列をコンパイルし、クラスから static なメソッドを呼んで値を返す方法を模索していました。

以下のプログラムと解説を参考にしています

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;
  }
}

2012-06-10

R Commander による Kaplan-Meier プロットの作図手順

R のインストール

Windows の場合
- Download R-3.2.3 for Windows. The R-project for statistical computing.
Mac の場合
- R for Mac OS X
- tcltk もインストールしておきます; R for Mac OS X - Development Tools and Libraries
- 参考情報
  - http://cse.niaes.affrc.go.jp/minaka/R/Rcmdr_for_MacOSX.html
  - Snow Leopardならクリックだけで美麗R Commander：Taglibro de H：So-netブログ
http://www.okada.jp.org/RWiki/?R%20%A4%CE%A5%A4%A5%F3%A5%B9%A5%C8%A1%BC%A5%EB に詳しい情報がありますので、困ったら参照しましょう

R Commander のインストール

install.packages("Rcmdr")

パッケージインストール画面の表示

ミラーサーバーを選択 (なるべく近そうなサーバを選ぶと良いでしょう)

インストールするパッケージの選択

インストールエラーが表示された場合は、もう一度同じ操作を行ってみてください。

R Commander プラグインのインストール

プラグイン

複数

install.packages("RcmdrPlugin.KMggplot2")

R Commander の起動

library("Rcmdr")

RcmdrPlugin.KMggplot2 0.2-0 以降では、R Commander の起動とプラグインのロードを同時に行うことができます。
以下のコマンドを実行するか、パッケージのロード画面から "RcmdrPlugin.KMggplot2" を選択します。

library("RcmdrPlugin.KMggplot2")

ロードが終了したら、5 番目のステップ (プラグインの読み込み) を実行する必要はありません。

パッケージのロード

初回起動時は、以下の様な画面が表示されるはずなので、両方共そのまま [はい] を選択します。選択すると必要なパッケージが自動的にインストールされます

Rcmdr の起動画面

R Commander プラグインの読み込み

プラグイン

プラグインのロード

再起動を促されるので再起動します

データの読み込み

イベント変数をイベントありの場合が 1、打ち切りの場合が 0 でコーディングした形式のデータを準備し、クリップボードにコピーしておきます

データ、データのインポート、テキストファイル・・・を選択します

データセット名を指定し、変数名はあり、クリップボード、タブをチェックし、OK をクリックします

作図

KMggplot2、Kaplan-Meier プロット... を選択します

時間変数、イベント変数、必要に応じてその他のオプションを設定し、OK をクリックします。グラフをファイルにセーブしたい場合はグラフの保存にチェックを入れておきます

実行結果です

2012-05-30

記号微分

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.)

deriv では基本的な演算子と、いくつかの単変数の関数ぐらいしか使えないため、被微分関数を作るときに注意が必要です。

2012-05-19

RcmdrPlugin.KMggplot2_0.1-0 is on CRAN now

R-bloggers

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.

`theme_simple()`

The theme_simple() is a simple black and white theme.
Top and right borders and both grids are removed from the theme_bw().

A list box for `facet_xx()` functions

A list box for facet_xx() functions was added.

Triad sou.

Recoding a factor variable in a data.frame by the levels() function

Recoding `z` as 1 = A, 2-6 = B, 7 = C

Result

Recoding `z` as 1 = "1", 2-6 = "2-6", 7 = "7"

Result

Penalized quasi-likelihood について

PQL

感想

参考文献

Log-gamma distribution

Note

Reference

Java Compiler API を使って何かしてみよう

ポイント

2012/06/25 追記

2012/06/26 リフレクションメモ

R Commander による Kaplan-Meier プロットの作図手順

記号微分

正規分布の対数尤度関数を記号微分してみよう

注意

RcmdrPlugin.KMggplot2_0.1-0 is on CRAN now

NEWS

Changes in version 0.1-0 (2012-05-18)

`theme_simple()`

A list box for `facet_xx()` functions

Confidence intervals for Kaplan-Meier plots

Violin plots

Stacked bar charts

Univariate plots at diagonal positions (`ggplot2::plotmatrix`)

Related Posts

Recoding z as 1 = A, 2-6 = B, 7 = C

Result

Recoding z as 1 = "1", 2-6 = "2-6", 7 = "7"

Result

PQL

感想

参考文献

Note

Reference

ポイント

2012/06/25 追記

2012/06/26 リフレクションメモ

正規分布の対数尤度関数を記号微分してみよう

注意

NEWS

Changes in version 0.1-0 (2012-05-18)

theme_simple()

A list box for facet_xx() functions

Confidence intervals for Kaplan-Meier plots

Violin plots

Stacked bar charts

Univariate plots at diagonal positions (ggplot2::plotmatrix)

Related Posts

Recoding `z` as 1 = A, 2-6 = B, 7 = C

Recoding `z` as 1 = "1", 2-6 = "2-6", 7 = "7"

`theme_simple()`

A list box for `facet_xx()` functions

Univariate plots at diagonal positions (`ggplot2::plotmatrix`)