VISA手続き備忘録

USのビザ手続きに関する備忘録。USのビザの手続きはすごく面倒です。必要書類はこちらにまとめてありますが、手順がわかりにくいので、こちらにまとめてみます。私の場合、海外学振のプログラムで(J1ビザ)2017年4月に渡航したときのタイムラインになります。

SPONSOR LINK

渡航前

2016年10月 DS-2019の手続きを開始する

まず受け入れ先の担当者と連絡をとり、DS-2019と呼ばれる書類(米国における長期滞在証)を手に入れます。DS2019の申請にあたり、米国滞在中のスポンサーからのレターが必要になりますので、それらが揃い次第コンタクトをとるのが良いかと思います。私の場合、2016年の10月末にスポンサーからレターをもらい、即座に連絡をとって準備を始めました。しかし、大学の国際担当部署が動かず、結局12月末にDS-2019が発送され1月上旬に届くという結果になりました。なお、2015年から英語能力に関わる証明が必要になっています(大学によっても基準が違うようです)。TOEFLやTOEICのようなものでスコアを持っていない場合、大学教員との面接をパスする必要があるので、早めの準備が必要です。

2017年1月中旬 DS-160の入力

DS-160は、渡航者の個人情報などを入力するフォームです。かなり長いストレスのたまる作業ですので、時間に余裕をもって取り組むことをお勧めします。最後のほうで*SEVIS IDと*Program Numberの入力があるのですが、これらはDS-2019がないと入力できません。そのため、DS2019を早めに手配するのが重要です。DS-2019が届かない場合、SEVIS IDは「N000000」はとし、Program Numberは受け入れ先の担当者に問い合わせて入手するという手段があるようです(参照)。

*SEVISとはStudent and Exchange Visitors Information Systemの略で、留学生(F-1:一般留学生、M-1:職業トレーニング生)や交換留学生(J-1)の監視・取締りを目的としたシステムのことです。その登録(支払いアリ)にはDS-2019に記載されている情報が必要。
*Program Numberとは、DS-2019に記載されている交換留学プログラムの識別番号になります。

2017年1月下旬 ビザ申請料金の支払いおよび面接予約

DS-160とDS-2019がそろえば、あとはビザの申請料金を支払い、面接予約(こちら)とることになります。ここまでくればだいぶ楽になりますが、この面接予約には少し注意が必要です。というのも、主要都市以外(東京など)の大使館では、それほど頻繁に面接が行われていません(私が面接した札幌では月2回程度)。そのため、そもそも日程を調整しずらいうえに、面接枠が埋まるのが早いです。東京なども面接枠が埋まりやすいことに変わりはないと思いますので、これまた早めの準備が肝心です。なお、申請先の大使館によると思いますが、面接前にパスポート以外の書類を事前郵送することになる場合があります。札幌大使館の場合、面接予約をすると空のレターパックが送られてきたので、そのなかに書類を入れて返送しました。

*提出書類の中に証明写真の提出が含まれるのですが、この規格に対するこだわりが半端ないです。すこしでも企画からはずれると(背景が白でない、など)自腹で撮り直しになるので、最初からビザ用の写真をとっている写真屋さんで撮ることをお勧めします。

2017年2月上旬 面接

札幌大使館の場合、面接は英語で行われていました。ですが、私の時の担当者の方は大変フレンドリーな方で、日本語もペラペラでした。アメリカ行くんだから英語の練習をしよう、くらいのノリで英語面接にしていました。何をしに行くんですか?とかその辺の普通のことを聞かれるので、研究分野名などに言及しておけば問題ないと思います。面接が終わると提出した書類一式は返却され、パスポートだけ後日郵送されてきます(ビザが添付される)。私のときは1週間強で届いたと思います。

渡航後

入国時

入国時には、「ビザ」と「DS-2019」の両方を入国管理官に見せる必要があります。このとき、キャリーケースなどにDS-2019を入れてしまうと別室行きになるので(最悪入国拒否?)、「ビザ」と「DS-2019」のどちらも携帯して飛行機にのるようにします。

入国後

J1ビザの場合、プログラム開始予定日から30日以内に米国に入国する必要があります。ただし、渡航者が大学の国際部署に入国の届け出を出し、大学から管理局に連絡がいくので、実際には10日~20日以内と考えておいた方がよさそうです。この期限に関しては受け入れ先の担当者に聞くのがよいと思います。ついたばかりのころは右も左もわからないので、余裕をもって到着しておくことをお勧めします。

入国後の住まい

住まい探しが最初の難関です。たいてい、大学のHPなどに居住先などに関するリンクが張ってあり、そこから探すのも一つの手です。私の場合最初の1週間は一時居住向きのアパートを探し、2週間ほどそこに滞在しました。1日30ドル前後で泊まれるので、この間に現地のアパートを見て回る、というのが現実的かと思います。私の場合、その後に教授があっせんしてくれたため、アパート探しは何も苦労せずにすみました。現地の院生とルームシェアから始まり、現在はスペイン人のポスドクとルームシェアしています。日本人ではないですが、人によっては最初の数週間はボスの家に転がり込み、その間にアパートを探している人もいるようです。

一時出国

一時出国の際は、DS-2019の右下にあるTravel validationという欄に国際部署の担当者によるサインが必要になります。これがないと再入国の際に大きな問題となり、最悪入国拒否される場合もあるようです。こちらに怖い経験談がありますので、一度目を通すとよいかもしれません。私も以前、このTravel validationの存在に出国前日に気づき、なくなく国際学会参加の全日程をキャンセルする、という悲しい事態に陥りました。学生にはオリエンテーションなどで広く周知されているようですが、ポスドクには明確に周知されない場合もあるようですのでお気をつけて。

濃淡のある信用区間のプロット

WinBUGSやJAGSを用いたベイズ推定を使うと、各パラメーターに対する大量のMCMCサンプルが得られます。このMCMCサンプルを使って95%信用区間を計算し、それをもとに議論しますが、こうした任意の信用区間ではなく、事後確率がどのような分布をしているのかのほうを見せたい場合もあります。そんな時に便利なパッケージdenstripを紹介します。

SPONSOR LINK

濃淡のあるパラメーター推定結果の図示

ここでは、行列XにパラメーターのMCMCサンプルの入っているとしてみましょう。各列が各パラメーター、各行がMCMCサンプルに対応します。

X <- matrix(NA,nrow=1000,ncol=3)
for(i in 1:ncol(X)){
  X[,i] <- rnorm(1000,rnorm(1,2,2),runif(1,0.1,1))  
}
colnames(X) <- c("p1","p2","p3")
> head(X)
           p1        p2       p3
[1,] 2.527981 0.9171227 4.432790
[2,] 3.007233 2.9618068 4.689230
[3,] 1.891604 2.8134809 4.616570
[4,] 2.741595 0.8363671 4.741110
[5,] 2.602177 3.3400819 4.502327
[6,] 1.603801 2.5635456 4.356336

このデータを使って、事後分布と比例するグラデーション付きのパラメーター推定値を図示してみます。手順としては、まず空のプロットを呼び出し、そこにdenstrip関数をつかってストライプを図示していくことになります。

library(denstrip) #load package
plot(0,type="n",xlim=range(X),ylim=c(1,3)) # make plot region
for(i in 1:3) denstrip(X[,i],at=i,mticks=median(X[,i]),mcol="white")
# at = i is y-axis values density stribe to be plotted
# mticks is the x-position for median or similar
# mcol is the tick color to be used

出来上がりはこんな感じです。濃淡は、特に指定しなければdensity関数から求められる確率密度カーネルを基に図示されます。

また、以下のようにすれば見栄えも調節できます。

plot(0,type="n",xlim=range(X),ylim=c(1,3), ann=F, axes=F)
for(i in 1:3){
  denstrip(X[,i],at=i,mticks=quantile(X[,i],c(0.025,0.5,0.975)),
            mcol=c("black","white","black"))
}#adding 95% credible intervals
box(bty="l")
axis(1); axis(2,at=c(1,2,3),labels=colnames(X),las=2)

濃淡のある線形回帰の信用区間

パラメーターの推定結果に加え、線形回帰の推定結果についても濃淡をつけて図示することが可能です。基本的な手順は一緒ですが、denstrip関数の代わりにdensregion関数を使います(densregion関数についての詳細はこちらを参照)。まず、説明変数xについて、線形モデルによる予測値のMCMCサンプルが入った行列を用意します。先ほどと同じように、列が予測値と対応しており、行がMCMCサンプルの繰り返しになります。

> head(MCMC[,1:3])
     ypred[1] ypred[2] ypred[3]…
[1,] 0.336092 0.333215 0.330362
[2,] 0.384028 0.379466 0.374958
[3,] 0.469216 0.467234 0.465260
[4,] 0.259841 0.256886 0.253964
[5,] 0.350182 0.347407 0.344654
[6,] 0.597118 0.593395 0.589696
…

このデータをもとに、下記のように整形していきます。

est <- apply(MCMC,2,median) # median estimate
x <- seq(min(P),max(P),length=100) # range of explanatory variable
y <- seq(0,1.6,length=512) # range of y-axis; probability density will be estimated in this region
z <- matrix(nrow=length(x),ncol=length(y)) # empty case for density values
for(i in 1:length(x)) z[i,] <- density(MCMC[,i],from=min(y),to=max(y))$y # insert density values

library(denstrip)
plot(0,type="n",ylim=range(y),xlim=range(x), ann=F, axes=F)
densregion(x,y,z, colmax="red", gamma=0.8,scale=0.75,nlevels=75)
lines(x,est,lwd=2,col=grey(0,0.65))
axis(1);axis(2,las=2);box(bty="l")

こんな感じになります。

魚釣りをする二枚貝

ときに生物は、人の想像も及ばないような進化を遂げていることがある。河川や湖沼にすむ淡水二枚貝のイシガイ類(Unionoida)は、その一つだと思う。

イシガイ類は世界で約800種が知られており(Strayer 2008)、そのほとんどが寄生生態をもつ。親貝自体はよく見る普通の貝なのだが、その幼生は魚(両生類の場合も)に寄生する。雌貝の中で成熟した幼生は水中へ放出され、周囲の魚に寄生する。寄生期間の長さは種によって様々だが、数週間から数か月であることが多い。特定の魚種しか利用しない貝(スペシャリスト)もいれば、選り好みしない貝(ジェネラリスト)もいる。なぜ寄生生態が進化したのかについては諸説あるが、移動分散のためなのはほぼ間違いないであろう。二枚貝は移動能力が低いので、魚にくっつくことで遠くへ移動できるという寸法である(Terui et al. 2014, 2015)。こうした寄生生態をもっているというだけでも十分驚きである。

しかし、下の写真、ちょっとしょぼい小魚に見えると思うのだが。。。

これをすこし引いて写真をとるとこうなる。

そう、この二枚貝は、まさに「ルアー」を作り出したのである。この二枚貝はLampsilis cardiumという種で、オオクチバスなどの肉食魚を宿主とする二枚貝である。ルアーの陰には寄生する準備のできた幼生がたんまりと仕込まれている。このルアーをエサと勘違いしてオオクチバスが攻撃した瞬間、”ぶわっ”と幼生を放り出して寄生するのだ。二枚貝とは思えない精巧な擬態と寄生戦略である。

川に潜ると、こんな感じでディスプレイしている様子を写真に収めることができる。いずれもL.cardiumである。

擬態は小魚に限られたものではない。ザリガニ、巻貝、イトミミズ、はてには「釣り糸まで作り出し、その先でルアーを泳がせる」二枚貝までいる。こうした寄生戦略の多様性はUnio Galleryにまとめられており、これを見るだけで壮大な進化の歴史に思いを馳せてしまう。いったいなにが、この二枚貝の異様な進化を引き起こしたのだろうか。いまだ謎に包まれたままである。

イシガイ類の多様性はミシシッピ川流域が群を抜いており、200種を超える(Haag 2012)。そして、極端な擬態が観察されるのも(今のところ)この地域だけだ。日本にも17種のイシガイ類がいるが、疑似餌を進化させた種はいない。そして、同じ北アメリカでも、ミシシッピ川流域を離れると、途端にその多様性は著しく低下し、疑似餌を持つ種はいなくなる(この辺は不正確;少ないことは確か)。ざっくりとした仮説としては、最終氷河期に氷河におおわれなかったため、多様性が高いとされている(これは二枚貝に限った話ではないけれども)。しかし、それだけではまったくロマンがない。なにかもっと面白いことが起きているはずである。そう思ってアメリカに来て、もう半年がたつ。いくつか面白いことはありそうなことが分かってきたが、それを実証するためのアプローチが思いつかない。苦悩は続く。

※残念ながら、イシガイ類は世界的に減少傾向にあり、北アメリカでは70%の種が絶滅危惧種あるいは要注意種になっている。日本でも17種のうち13種が絶滅危惧種もしくは準絶滅危惧種である。

Strayer 2008 Freshwater mussel ecology: a multifactor approach to distribution and abundance. University of California Press

Haag 2012 North American Freshwater Mussels: Natural History, Ecology, and Conservation. Cambridge University Press

Terui et al. 2014 Dispersal of larvae of Margaritifera laevis by its host fish. Freshwater Science 33:112-123

Terui et al. 2015 A “parasite-tag” approach reveals long-distance dispersal of the riverine mussel Margaritifera laevis by its host fish. Hydrobiologia 760: 189-196

Ecology主要ジャーナルの論文の長さ

長い論文がよい?それとも短い論文がよい?

内容はともかく、やたら長い論文を読むとものすごい時間がかかるので嫌気がさすことがある。そんな話をスペイン人のルームメイトに話をしたら、彼はものすごい長い論文を書くタイプの人だった(げっ)。はてさて、私は結構短い論文を書くタイプの人間だと思うのだけど、それは実際のところどうなのだろうか、とちょっと気になった。そこで、Ecologyの主要ジャーナルについて論文の長さを調べてみた。

IntroductionとDiscussionの比較

とりあえずIntroductionとDiscussionの比較からしてみよう。感覚的にはDiscussionのほうがやや長くなりがちだと思うのだが、いったいそうなのであろうか。以下の7誌(Ecology, Ecology Letters, Ecosystems, Journal of Animal Ecology, Oikos, Oecologia, Proceedings B)から直近の10報を無作為に取り出し、それらの各セクションの長さの統計を取ってみる。

この図のx軸は、Introduction(黒線)とDiscussion(赤線)のワード数を表している。Journalに関係なく集計したもので、縦軸の高いところほどそのワード数(x軸)の論文の頻度が高いことを意味する。こうしてみると、やはりIntroductionよりもDiscussionのほうが概して長い傾向にあることがわかる。それぞれの中央値は、Introductionが876 words、Discussionが1454 wordsであった。これらの値を超える場合、それは全体の50%分位点より高いことになるので、比較的長い論文ということになるだろう。ちなみにIntroductionに限定すれば、1000 wordsを超える論文はほとんどなかった。

Journal別にみてみる

さて、文章の長さには、Journalによる傾向はあるのだろうか。まずはIntroductionからみてみよう。

上の図は、白丸は平均値、バーがSDを表している。Journal of Animal EcologyからProceedings Bの順に短くなっていた。Royal Society系は短い論文が多く、BES系はぐだぐだ長い論文が多い印象を受けていたので、この結果は納得がいく。

次はDiscussionについてみてみる。

EcosystemsからOecologiaの順に短くなっていた。Ecosystemsは極端に長く、Oecologiaがやや短くなっていた。他はおおむね1500 wordsのあたりに落ち着いていた。1500 wordsが適度なのだろうか。私は長すぎると思うのだが。。。

自分の論文と比較する

自分が書いた直近2論文と比較してみる。
Ecosystems – Intro: 694 , Discussion: 909
Proceedings B – Intro: 772, Discussion: 1011

いずれも中央値より全然短い。やはり短く書くタイプだった。

意味するところ

上の統計の意味するところはそれほど深くはないだろうけれども、論文を書く際のちょっとした参考値にはなるのではなかろうか。すくなくとも上記のJournalに投稿され、受理された論文については、上に述べたword数くらいの説明は必要とされることが多いということだろう。あまりにもこれらの統計値から外れる場合、極端に説明不足か、過剰な説明になっている可能性がある。ただし、各Journalから10報ずつしかサンプルを取ってきていないので、サンプル不足による統計力不足は否めないのであしからず。

読みやすい英語を書くための四つの”べし”

英語しんどい

論文を書くにあたって、日本人がぶち当たる最初で最後の壁は「英語」である。これはどうしようもない壁であるが、どうにかしないといけない壁でもある。そこで、大学院時にネイティブの先生に教えてもらった読みやすい英文を書くための覚書(+自身の経験則)を(自戒を込めて)残しておこうと思う。

1.主語と動詞のつながりをすっきりと見せるべし

ついついやってしまいがちなのが、なが~い主語。これは全然だめらしい。私がよく指摘されたミス(意味は通るけどわかりにくい)は、クラウズ(whichなどの関係代名詞)をつかって主語を長々と修飾してしまうことだった。なぜダメなのか。その文の主役である「主語と動詞」が初見でみつけづらいからだ。でもそれだと正確な英文書けないよ!という場合。それはそもそも適切な文構造を選び損ねている可能性が高い。能動態・受動態のどちらがよいのか、”It is”から始まる文にした方がよいのか、同じ意味の動詞でも自動詞と他動詞をスイッチすることでどうにかならないか、、、etc。「主語-動詞」のつながりがパッとみえる文構造をひねり出す。これだけで、「論文の可読性」はグンと上がる。たぶん。

2.クラウズはパラグラフ当たり1個くらいに絞るべし

1とも関連するが、クラウズの多用は厳禁である。日本語と異なり、クラウズを使うことで補足情報をどんどん付け足せるのが英語の特徴的な部分の一つである(と私は思っている)。だけどそれに頼ってはいけない。なぜか。クラウズのあまりにも多い英文は、いったい真に重要な情報は何なのかわかりにくいからだ。たぶん、クラウズを使いまくっているうちに、書いてる本人も何が大事か忘れ去っている。

3.長い一文は避けるべし

長い一文、とくに意味のない重文は絶対に避けるべきだ。むやみにAndやAlthoughのような接続詞を使って重文にすると、それだけで「ん?」と思われる可能性が高まる。重文は、重ねることでさらなる理解が促される場合にのみ使う。Becauseなどの理由を説明する重文はOKな場合が多い。

4.同じ表現の繰り返しは避けるべし

これは上記の3つと比べると重要度は低いが、できるだけ同じ構文やワードを連続して使うのは避けたほうが良い。これは言語的な特徴かもしれないが、英語は極端に同じ表現の繰り返しを嫌う傾向にある(気がする)。これを達成するためには、とにかく同じ意味でも複数のワードや構文を頭にストックしておく必要がある。これは今すぐどうにかなる問題ではなく、普段からどれだけ英語論文を(意識して)読んでいるか、ということに終始する。

だがしかし、冠詞はいつまでたってもわからん。

ジャックをジャァァックなどと呼んではいけない

アメリカにきてそろそろ半年がたつ。やはり一年目は苦労するだろうなぁなどと思っていたが、まぁ英語の苦労もそこそこに、一年目からデータが取れたので少し安心した。最悪の場合も想定し、データが全く取れずとも何とかなるバックアッププランも用意していたのだが、そこまでは至らなかった。よかった。

名前の発音

しかし意外なところで苦労した経験がある。それは名前の発音である。現在私が所属しているラボのボスは、ジャック(Jacques)という方である。EEB(Department of Ecology, Evolution and Behavior;現所属)のコーヒータイムなどで「ここではだれと仕事しているの?」と聞かれるので、ジャックだよ、と当然答える。が、十中八九は「え、だれ?」といわれてしまう。最初は、あれ、部署も同じだし知らないはずがないよなぁ、と首をかしげていたのだが、しばらくたって分かったのは、私の発音が「ジャァァック」のように変に間延びした発音になっていたためらしい。ジャックはジャッッ”クッ”くらいの心持で発音しないといけない。

う~ん、英語って、むつかしい。

不運なRejectを避けるために

主著の論文が10本を越えてから、やっとIF > 4以上の中堅雑誌(Ecologyの分野では)に載せることができるようになってきた。それまでもレビューに回ることはあったが、一度とてアクセプトにこぎつけることはできなかった。さて、その勝敗の差を分けたのはなんなのであろうか、と自分なりに整理してみる。ただし、以下は十分良い結果が得られている前提で、いかにその結果を美味しく料理するか、という段階であるのであしからず。

一般性の担保

これは言わずもがなかもしれないが、IF > 4になると、特定の生態系、あるいは特定の分類群にしか当てはまらないことが明らかな現象は、Editor Rejectを食らう運命にある。ただ、普通は一種、よくても数種しか扱えないと思うので、いかに現象論として一般化できるのかをイントロとディスカッションで大きく、かつ説得力のある形で打ち出せるかが勝負どころ。大きく出るだけであれば簡単だが、やはり説得力を持たせるためにはかなり綿密な既存文献のレビューが必要になる。

新奇性の見せ方

新奇性があるにも関わらず、ぐだぐだと長い文章を書き連ね、新奇性がイントロの後のほうに出てくるのはよくない気がする。重要な知見の欠落を早い段階で指摘し(できれば2段落目、遅くとも3段落目)、そのギャップをいかにインパクトのある形で見せれるかが違いを生み出す気がする。

手法の説明

ここは見落としていたポイントだった。適切な手法を使っていようと、それが読者に伝わらなければ、それはトンチンカンな理由からリジェクトを食らう原因となる。私はベイズモデルをよく利用するのだが、その利点をあまり説明しなかったがために、「古典的なモデルを使ったほうが良い」という理由からリジェクトされた経験が何度かある。最近は、GLMMなどと比較して、こうした大きな利点がある、というのを明記するようにしており、この理由によるリジェクトはグッと減った。何度も「おれの方が正しいのに!」と思ったが、思い返してみれば「自分の手法の正しさをわかりやすく記述していない」ことが原因にあったのだ。自業自得であった。

結果の書き方

これも盲点だったのだが、単に「~の有意な効果が得られた」ではいけない。「コントロールと比べ、3~4倍もの~が生じることがわかった」のような、定量的なインパクトを示す記述が必要。いわゆる普通のジャーナルとトップジャーナルのResultの書き方を見比べると、この点は雲泥の差があることがわかる。印象が全然違うのだ、きっと。

弱点のカバーの仕方

どんな研究にも弱点はある。だが、弱点をただ認めるだけでいいのは下のレベルのジャーナルだけのような気がする。いかに弱点を「弱点らしからぬフリ」で守り抜くかが実は大事なテクニックであるように感じる。いかにさりげなく、だけどしっかりと研究の穴をガードするのか、これはなかなか難しく、数をこなさないとできないのだろう(まだ検討中)。

Supporting Informationの有効活用

想定されるネガティブなコメントで対処できるものであれば、面倒くさがらずにSIのなかで「こういう場合も想定されるが、それでも同様の結果が得られる」といったような記述を加えておく。上のJournalになるほど、ちょっとしたネガティブなコメントが入るだけで、一発リジェクトの危険性が高い。そうなる前に、予防線を張っておくのが無難である。

英語力

とにかく文献を読み、片っ端からいい表現を盗み取るに限る(文自体コピペはいかんけど、あくまで英語表現)。つたない英語表現より、読みやすい表現のほうが絶対アクセプトされやすい。また、経験上、つたない文章は、英文校閲による改善も限界がある。一定ライン(ゼロベンでTOEIC700~800点くらい?)を越えた英文は、校閲者がきちんと意図をくみ取ってくれるので、よりよい英語表現を的確にSuggestしてくれる(これは自分の英語力の改善と共にひしひしと感じている)。結局自分自身の英語能力を上げるのは必要不可欠で、常日頃から意識して鍛えるべき能力という他ない。

行列操作

Rの行列操作になれる

Rの行列指定になれると、いろいろと作業が簡略化できます。備忘録もかねて、整理したいと思います。

数字による行列指定

Rでは「行列[行番号,列番号]」の形で行列が指定されます。例として、3×3の行列を作り、特定のデータを行列の番号で指定してみます(Xは行列ではなく、データフレームでも可)。

> X <- matrix(rnorm(9,0,1),3,3)
> X
          [,1]       [,2]       [,3]
[1,] 2.3587491 -0.8524698  1.8697275
[2,] 0.8162371  0.4589377 -0.3638491
[3,] 1.8183986  1.2378938 -0.7109136

> X[2,3]
[1] -0.3638491

複数のデータを同時に指定する場合は、”1:3″(1から3まで)や”c(1,3)”(1と3を指定)などを使います。また、空欄”[1,]”とすると、その空欄の部分(行もしくは列)のすべての要素が指定されます。

> X[c(1,2),1]
[1] 2.3587491 0.8162371

> X[1:3,1]
[1] 2.3587491 0.8162371 1.8183986

> X[1,c(1,2)]
[1]  2.3587491 -0.8524698

> X[1,]
[1]  2.3587491 -0.8524698  1.8697275

特定の行や列を取り除きたいときは、数字を負にします。

> X[,-2]
          [,1]       [,2]
[1,] 2.3587491  1.8697275
[2,] 0.8162371 -0.3638491
[3,] 1.8183986 -0.7109136

> X[,c(-2,-3)]
[1] 2.3587491 0.8162371 1.8183986

要素名による行列指定

要素名を指定することも可能です。

> colnames(X)<-c("a","b","c")
> X
             a          b          c
[1,] 2.3587491 -0.8524698  1.8697275
[2,] 0.8162371  0.4589377 -0.3638491
[3,] 1.8183986  1.2378938 -0.7109136

> X[,"a"]
[1] 2.3587491 0.8162371 1.8183986

複数ある行のうち、特定の列名(行名)に一致する列番号を取り出したい場合は、colnames(rownames)を使います。複数の列名(行名)に一致する番号を知りたい場合は、”==”ではなく”%in%”を使います。

> colnames(X)
[1] "a" "b" "c"

> which(colnames(X)=="a")
[1] 1

> which(colnames(X) %in% c("a","b"))
[1] 1 2

> X[,which(colnames(X) %in% c("a","b"))]
             a          b
[1,] 2.3587491 -0.8524698
[2,] 0.8162371  0.4589377
[3,] 1.8183986  1.2378938

> X[,-which(colnames(X) %in% c("a","b"))]
[1]  1.8697275 -0.3638491 -0.7109136

複数ある行のうち、特定の列名(行名)に一致する列番号だけをのぞきたい場合は、”!=”も使えます。ただし、これは指定する行もしくは列が一つだけの場合に限られるので、「-which(colnames(X) %in% c(“a”,”b”))」などとするほうが汎用性はあるように思います。

function関数の使い方

function関数

Rで複数の関数を組み合わせて値を算出しなければならないとき、function関数で定義してしまうと楽です。しかし、いったい何をやってるのかわかりにくいところもあります。ここでは、function関数を使った関数の定義について簡単に紹介します。

変動係数

変動係数(CV)とは、標準偏差を平均で除したもので、変動性の比較などでよく使われる指標です。RのデフォルトでCVは定義されていないのですが、比較的よく使うものなので、例として挙げてみます。まず、function関数の仕組みですが、引数の指定と、その引数を使った計算式の指定、という形になっています。

CV <- function(x){ sd(x)/mean(x)}

#example
> x <- rnorm(100,10,10)
> CV(x)
[1] 1.042615

上記のfunction()のかっこ内に引数(今回はx)として何を使うのか指定し、{}内に引数を使った計算式を定義します。

ベクター中の0要素をカウントする関数

少し込み入ったものとして、ベクター中に0がいくつかあるのかカウントし、そして0が何番目の要素かを返す関数を定義してみます。

zerocount <- function(x){
 n_zero <- length(which(x==0))
 id_zero <- which(x==0)
 return(list(n_zero=n_zero, id_zero=id_zero))
}

#example
>x <- rpois(100, 2)
> zerocount(x)
$n_zero
[1] 17

$id_zero
 [1] 12 16 21 22 25 31 35 40 43 67 74 77 78 79 86 87 95
 

こちらでは、返す要素が二つあるので、返り値をリスト化して返すよう指定しています(4行目)。

応用

発展的には、乱数発生などのプロセスを組み合わせることで、簡単なシミュレーションモデルを作ることができます。そのほか、自作の作図関数なども作ることができます。興味のある人はぜひトライしてみてください。

最頻値の推定関数

最頻値

平均、中央値、分散など、確率変数を特徴づける指標はいくつもあるけれども、多くのものについてはRで標準実装されている。しかし、地味に知りたい「最頻値」なるものは実装されていない。

最頻値の推定の仕方

離散値の場合は以下のコードで簡単に推定できる。

x <- rpois(100,10)
names(which.max(table(x)))

しかし、サンプルサイズが比較的小さい場合や、確率変数が連続値の場合にはとてもいい方法とは思えない。そこで、table関数の代わりに、density関数をつかった関数を定義する。

mfv <- function(x,from=min(x,na.rm=T),to=max(x,na.rm = T)){
  if(is.numeric(x)==F)stop("use numeric")
  tmp <- density(x,from=from,to=to)
  Mode <- tmp$x[which.max(tmp$y)]
  return(Mode)
}

# example
> x <- rnorm(1000,10,20)
> mfv(x)
[1] 12.16744

density関数は、確率変数の確率密度(のようなもの;実際はカーネル)を変数の傾度にそって推定してくれるもの。上で定義したものは、このカーネルの値が最も高くなる確率変数xの値を返す、という関数になる。細かいけど、いちいち計算するのは面倒なので、関数化してしまうのが楽。。。ちなみに、fromとtoに値を指定すると、その範囲の中だけでカーネルを推定する。例えば、カウントのような正の値しかとらない変数の最頻値を知りたいときは、from=0とすることで、正の実数の範囲の中だけで最頻値を計算してくれる。また、確率など0-1しかとらない、とわかっているのであれば、from=0、to=1とする。

だが、そもそもdensity関数がなにをどうやって推定しているのかようわかっていないので、正確ではない可能性がある。確認次第追記します。