（cf. Rいろはの入り口 http://nhkuma.blogspot.jp/p/r_5.html)
# 注：作成者は統計の専門家ではありません。内容には十分な注意を払っているつもりですが、限界があることを理解した上で参照してください。感想・苦情・間違いのご指摘、歓迎します。
（注：使用環境によっては、スペースやクオーテーションマーク"などが化けてコードがエラーになるかもしれません）
###### R基本操作 ######################
# 今日はRらしく、コードで行ってみます。１行１行、数式を書いて実行して進めていくソフトです。コード打つのを怖がらないでくださいね…
# Excel: 大型ショッピングモールの広大な平地の駐車場に停められた車を扱うような感じ
# （配置を覚えておかなければ迷子になる。何をするにも毎回停めてある場所まで行って探し出さなければ始まらない。方向音痴泣かせ…哀）
# R: 立体駐車場の入り口でIDを告げて預けられた車を呼び出すような感じ
#（こちらでIDさえ指定しておけば、配置や順番は問題にならない。ただし、現実の駐車場のような出し入れの煩わしさは無く、一瞬で預けられ呼び出せる！）
# コード打ちは試行錯誤が付きもの。なので、エディタ（ワープロソフト）に書き出し、保存しておくと利用しやすい。
# Rのこの画面に直接打ち込むか、コピペするなどして実行していきます↓
###### 基本の"き" ##############################
# シャープ以降はその行は無視される（メモ書きに便利）
# 以下、これ通りに打って、出力を確認してみよう（メモ書き部分は打たなくていいです^^;）
2 + 1 # 足し算
3 - 1 # 引き算
2 * 5 # 掛け算
9 / 3 # 割り算
# [1] というのは、単に１番目の出力結果であることを示す。そのうち、複数結果の例が出てきます。
# + 等の前後の半角スペースは付けなくてもいいが、付けた方が分かりやすい
# ただし、全角スペースは厳禁
a <- 3 # 変数を定義する（右辺で左辺を定義する）
# 矢印は < と - の組み合わせ。矢印の前後も半角スペースをいれるべし
a # 定義した後でaと打つと3を返す、以下同様に、定義したものを出力してみよう
majimun <- 4 # 変数名は自由（数字から始めてはいけない、スペースやハイフンは使用不可）
b <- "A" # 文字は""で囲む
d <- "haisai" # こんなのもOK
d2 <- "夜露死苦" # 全角も行けないことはないが、半角英数だけにしておいた方が無難
d3 <- "(-_-#)/" # こら！そこ、脱線しないように！
e <- c(1, 2, 3, 4, 5) # cで囲むと数字列（ベクトルという。"c"はcombineの頭文字を意味する）
f <- c(1:5) # 上と同じことを表す、コロンで結んでもよい
f + 1 # 全部に1を足す
f <- f + 1 # 矢印はイコールではない。同じ変数を使ってアップデート（再定義）できる。
f + c(5:1) # ベクトル同士の足し算
c(2:4) * c(3:5) # ベクトル同士の掛け算
rep(1, 10) # 1を10個
rep(c(1:5), 3) # 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5
rep(c(1:5), each=3) # 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5
seq(0, 10, by=2) # 0 2 4 6 8 10。0~10を2ずつ増やす
seq(0, 10, length=3) # 0 5 10。0~10を3等分
1e2 # = 1.0 x 10^2を示している
2.5e-10 # = 2.5 x 10^(-10)）
##### 文字列、文字列の操作 #########################
letters # 小文字のアルファベット一覧（下記の”データの操作”の要領で一部抜き出し可能）
LETTERS # 大文字のアルファベット一覧
month.abb # 月の3文字表記の一覧
month.name # 月のフル表記の一覧
paste(c("A", "B"), 1, sep=".") # "A.1"、"B.1"。文字の合成（sep=の中身はAと1の間に挟む区切り文字）
paste(letters[1:3], c(1:3), sep="", collapse=" + ") # "a1 + b2 + c3"。sepとcollapseは、区切り方の違い
substr("abcdef", 2, 4) # 文字の抜き出し（この例では2~4文字目を抜き出し、bcdとなる）
##### 行列（二次元のデータの表し方）################
matrix(0, 5, 5) # 行列（0を5行5列で）
matrix(c(1:16), 4, 4) # 1~16の4×4行列
v <- matrix(c(1:16), 4, 4, byrow=T) # 行方向に並べる場合（by row equal True）
colnames(v) <- c("sp1", "sp2", "sp3", "sp4") # ラベルを付ける
# cf. 行に名前を付けるのはrownames()
v # ラベルを付けたあとはこういう表示になる（実行例）
###### 基本的な関数 ##############################
# 平均 mean, 最大 max, 最小 min, 範囲 range, 標準偏差 sd, 分散 var, 大きさの順位 rank,
# データの個数 length, 行列の行数 nrow, 行列の列数 ncol
# 絶対値 abs, 指数 exp, 自然対数 log,
# 三角関数 sin, cos, tan # 円周率 pi を使用するのが基本：例 cos(pi/4) = cos(pi*45/180) である
# 正負の符号 sign, 小数点の丸め round, 小数点の切り捨て trunc
# 有効数字を指定して出力（, 桁数の数字を入れて使用） signif
# %/% 商（例：5 %/% 2）, %% 余り（例：5 %% 2; cf.数式で表す時：mod 2）
# 注：掛け算が入っている場合、直前の数字にしか適用されないので注意（両方に適用させたい場合：(y*5)%%2のように括弧で括る必要がある）
mean(f) # 平均値の算出の例
log(f + 1) # 対数変換をするなら（0データ混じりの場合、1を足しておけばエラーにならない）
sqrt(f + 0.5) # べき乗変換（0データ混じりの場合、0.5を足しておけばエラーにならない）
### 自作関数の利用 #############
# 標準誤差（よく用いる平均値の標準誤差）の関数はないので、それをse()関数として自作してみる
se <- function(x) sd(x) / sqrt(length(x)) # function(x)というのが、以降xの関数ですよという意味
# この場合、"x" 部分は他のアルファベットや単語でもOK。半角スペースは一応、空けておきましょう
se(f) # 使用例
# 関数の説明を見たい時
help(関数名) # ただし英語。必ずしも懇切丁寧ではないが…
?関数名 # help(関数名)と同様
??パッケージ名 # パッケージの概要と内包する関数一覧
example(関数名) # とやると、実行例が見られる。必ずしも分かりやすくはないが…
##### データの操作、範囲指定 ################################
# 上で定義したベクトル、
f <- c(2:6) # を例に（順に実行していればこの値になっているはず）。
f[2] # ベクトルの一部を取り出し（2番目を取り出し）
f[1:4] # ベクトルの一部を取り出し（1~4番目を取り出し）
f[c(1,3)] # ベクトルの一部を取り出し（飛び飛びで、1,3番目を取り出し）
f[f != 3] # ベクトルの一部を取り出し（値が3のもの以外）
f[-1] # ベクトルの一部を取り出し（1番目以外）
f[f >= 3] # ベクトルの一部を取り出し（3以上）
# cf. =と組み合わせる場合、不等号や!の位置は = の左側
# 上で定義した行列、
v <- matrix(c(1:16), 4, 4, byrow=T) を例に。
v[2,] # 行列の行を取り出し。[行番号, 列番号]を意味する
v[,3] # 行列の列を取り出し
v[2,2] # 行列の 1 点を取り出し
v[1:2, ] # 複数行を取り出し
rbind(v[2,], v[2,3]) # データの追加：上下で連結する（cf. 左右で連結は cbind。r: row, c: column）
which(v < 4, T) # 条件を満たす値の行（row）列（column）の番号（Tを除くとベクトル位置で返ってくる）
which.min(v) # ベクトルの最小値の位置（番号）を求める（cf. which.max # 最大値）
any(), all() # ベクトル要素の１つでも / すべてがカッコ内の条件を満たす時…if()との組み合わせが多そう？
##### 作業場所（データ階層）の指定（作業ディレクトリの指定）#########
# Rに読み込みたいデータファイルを置いた場所を作業ディレクトリに指定せよ
# Mac: メニューの"その他"→"作業ディレクトリの変更"。 command + d, command + d, returnで"デスクトップを指定できる。
# Win: メニューの"ファイル"→"ディレクトリの変更"
# cf. 作業ディレクトリの変更 http://cse.naro.affrc.go.jp/takezawa/r-tips/r/06.html
getwd() # 現在のディレクトリを確認する。この場合（）内は空白
setwd("でぃれくとり名") # ディレクトリの指定をする
##### データの準備・読み込み ##############################
## （cf. R解析用データの作り方）
# （予め用意したデータのファイルを読み込む場合）
# Excel などで csv 形式で保存（「カンマ区切りファイル」と出ているはず。Comma Separated Value のこと）
# 一行目を見出しとしたデータベース型を取ること
# データを "test.csv" というファイル名に保存し、作業ディレクトリに置き、そのデータを "test" としてRに読み込む場合のコード：
test <- read.csv("test.csv") # read.csv: データの読み込み
# あるいは、
test <- read.csv(file.choose()) # こちらはファイルメニューから選択（作業ディレクトリ外からも読み込み可能、こっちのほうが便利そうだ）
test <- read.csv(file.choose(), fileEncoding="cp932") # 日本語データを含むデータを読ませようとして "〜に不正なマルチバイト文字があります" が出た場合、このようにエンコーディング指定するとけっこう読めたりします
# test.csv のサンプルデータを用意しておきました（この出し方は理解する必要なし）
X <- rep(c(1:5),each=40) # rep(何を,何回) 繰り返す
F1 <- rep(rep(c("A","B"),each=20),5) # 文字を繰り返す場合は""で挟む
F2 <- rep(rep(c("A1","A2","B1","B2"),each=10),5)
Y <- rpois(rep(1,200), # rpois(いくつ,平均値) のポアソン乱数（整数値）を発生
rep(seq(20,200,by=40),each=40)*rep(rep(c(1,1.4),each=20),5)*rep(rep(c(1.0,1.2,1.4,1.6),each=10),5))
test <- data.frame(X,F1,F2,Y) # ひとつのデータセットに統合
# 一見、行列に見えるが、文字列と数値列が混在していても、それぞれ文字と数値として認識してくれる
# それがデータフレーム。Rでの解析の基本です。read.csvかdata.frame関数で作りましょう。
head(test, N) # データ内容の確認：上から N 行を表示（大量データの確認用）
# cf. tail() # データを下から表示
##### データの表記ルール ##################################
# 空白やスペースのみのセルは極力避ける（たいてい平気ですがエラーの元になりやすい）。ゼロなら0、欠損ならNAと入力する（NA：欠損、no data）
# 半角英数で記述すること（全角は極力避ける）。ピリオドやアンダーバーはいいがハイフンはダメ。
# アルファベット順に順序づけられる
#（e.g. January, Februaryと入れると、Februaryの方が先 → 頭に数字を付けると簡単。1January, 2February）
# 下記の「種の集計」のところにあるfactor関数を利用しても順序づけることが可能
# データ表記方法の詳細はこちらを参照
##### 欠損値（NA）の扱い方 ##################################
# Rでは欠損値をNAと表記しておく
#（
# しかし、関数にもよるが、NAが含まれたままではトラブルのもと。
# 例えば、mean(c(1,NA)), max(c(2,NA)) などはいずれも計算結果はNAになる。
# NAの取り扱いは案外厄介…下記のようにやるとよいでしょう
# cf. "!"は"それ以外"を意味します。
test[!complete.cases(test),] # 欠損値のある行を表示（箇所を見つけて手動で元データを修正する場合に使用）
for (i in 1:(ncol(test))) test[,i][is.na(test[,i])] <- 0 # 欠損値を0に置き換える（構造が複雑だが、”test”部分の書き換えだけで二次利用可能）
test2 <- na.omit(test) # 欠損値を含む行を除く
v[!is.na(v)] # ベクトルならば、このようにしてNAを除ける。
# meanやmaxなどでは、na.rm=Tという項で対処も可能だが、maxとminでNAのみのデータでは-InfやInfになる。!is.naを使っても0になる。NAの計算はNAになってほしいのだが…
##### パッケージのインストール #############################
install.packages("Rcmdr", dep=T) # パッケージをインストール（この場合はRコマンダーの例）、リストアップしてメモしておくとRを再インストールしたとき便利（定期的に最新にすべし）
update.packages("Rcmdr") # パッケージのアップデート
##### データの取り出し #######################
write.csv(test, "output.csv") # ファイルへ書き出し（データ構造により書き出せない場合あり）
# cf. Rの出力画面からコピペしたデータを加工するのに便利なサイト
# （単にコピペしただけだと、空白はスペースになってしまう））
# http://filter.webings.net/
# "連続半角スペース（1文字含む）→タブに変換" を選択、実行
# 最終手段はワープロアプリにコピペして、連続半角スペースを ,（カンマ）に変換し、拡張子をcsvに書き換えるとExcel上で復活させられる（多少ずれるでしょうけれど）
###### データ（データフレーム）の集計、ちょっとした計算など ##########################
test$F1 # testのF1だけを取り出し
test[test$F1=="A", ] # F1がAであるデータを取り出し
test[test$F1=="A" & test$F2=="A1", ] # F1がAで、F2がA1であるデータを取り出し
test[is.element(test$F1, c("A", "B")), ] # F1がAまたはBのデータを取り出し
unique(test) # 全く同じ行が複数ある場合、１データ行のみを残す（データ重複を無くす）
test[!duplicated(test$F1),] # F1の重複がないデータを残す（uniqueに近い用途）
# [test$F1, ] のようにいちいちtest$を付けるのが面倒な場合、いったんattach(test)とやると、単に[F1,]のように省略できるようになります。ただし後々、より複雑な処理をするにつれattachの使用はエラーの元になります。例えばGLMMのlmer関数ではattachは通用しません。
こういうときはwith関数をお勧めしています。with(data, lmer(~ formula))のように使えます。
###### ある列を基準にデータの並べ替えをする ######################
# Excelには対応関数無し（メニューから実行するしかない）
test[order(test$F2, test$F1),] # F2, F1の優先順位で全体行を並べ替え
test[order(test$F2, decreasing=T),] # 降順で並べ替え
###### ヒストグラム #####################
# 頻度別に集計しなくていい（Excelではアドインを使った大がかりな作業…）
hist(test$Y)
hist(test$Y, breaks=c(0:5)*100) # 区切りを指定することもできる
###### 種の集計に便利（必須？）#####################
# （Excelではcountifを使った繁雑な作業が必要…）
spp1 <- rep(c("sp1","sp2","sp3","sp4"),c(5,4,3,0)) # サンプルデータ
spp10 <- rep(c("sp1","sp2","sp3","sp4"),c(0,2,3,1)) # サンプルデータ2
table(spp1) # 種の集計：これだけでも行けるが、いない種はインデックスが作られない
species <- c("sp1","sp2","sp3","sp4") # これを解消するにはインデックスを予め用意する
spp2 <- factor(spp1, levels=species) # インデックスと照合（levelsの順で認識する）
table(spp2) # いない種についてはゼロとなる
unique(spp1) # spp1の要素をリストアップ（重複無しのリスト）
union(spp1, spp10) # 和集合。spp1、またはspp10に含まれる全種リスト
intersect(spp1, spp10) # 積集合。spp1、かつspp10に共通な種のみのリスト
setdiff(spp1, spp10) # 差集合。spp1にのみ含まれる種のリスト
combn(spp1, 2) # 組み合わせ。spp1から2種を取り出す場合（各列が組み合わせ）
# 各種の多様度指数、多様度の解析
# （以下、community1, 2というデータを仮定し、計算してみる）
community1 <- sample(1:5, 10, replace=T) # 10種の個体数を想定（1~5個体の範囲で）
community2 <- t(sapply(1:5, function(s) sample(1:5, 10, replace=T))) # さらに、行列版
library(vegan) # 生物群集解析のためのパッケージ
diversity(community1) # Shannon-Wiener index
# cf. 常用対数にする場合は, diversity(community1, base=2)
diversity(community1, index="simpson") # Simpson index (D)
diversity(community1, index="invsimpson") # Simpson index (1/D)
sim <- vegdist(community2, method="bray")
# 各行の群集間類似度を対戦表で表示（Bray-Curtis dissimilarity）
# 類似度（Bray-Curtis similarity）にするには全体を 1 から引く（1 - sim）
# cf. method="jaccard"：Jaccard類似度（在・不在データやランクの場合）
cruster <- hclust(sim, method="average") # 群平均法を用いたクラスター分析
plot(cruster, hang=-1) # custerを図示、hang=-1 で葉を揃える
# その他、MDS（多次元尺度構成法）なども"vegan"に含まれる
# vegan全般： http://cc.oulu.fi/~jarioksa/opetus/metodi/vegantutor.pdf
# veganで序列化法（MDS, CCA, RDAなど）：http://cran.r-project.org/web/packages/vegan/vignettes/intro-vegan.pdf
##### 集計やまとめて統計量を算出するのに便利な技 ###############
# 関数、の部分は、mean, sd, rankなど、なんでも
apply(test[,c(1,4)], 2, 関数) # 1:各行、2:各列に対し、関数を適用（この例では各列のmean）
aggregate(data=test, Y ~ F1, FUN=関数) # F1の各要素について関数を適用（excelのsumif等の拡張版）
aggregate(data=test, Y ~ F1 +F2, FUN=関数) # F1, F2の各要素について（無い部分はNAになる）
aggregate(data=test, Y ~ F1, FUN=function(x) sqrt(var(x) / length(x))) # SE を定義しつつ利用
sapply(1:5, exp) # 各値を関数に順に適用
sapply(1:5, function(x) x^2 + 2*x + 1)
# 各値を、x^2 + 2*x + 1 みたいに回帰式に当てはめるなど
combn(4, 2, 関数) # 組み合わせのそれぞれに対し関数を適用
##### 条件分岐＝プログラミングの初歩 #############
# 条件を満たすか否かで計算式を変える
# 上記で定義した
a <- 3 # を用いて…
if (a > 0) { "positive" } else { "negative or zero" } # 0より大か否かを判別
if (a > 0) { "positive" } else { "negative or zero" } # 0より大か否かを判別
if (a > 0 & a < 2) { "0 < a < 2" } else { "no" } # 0 < a < 2か否かを判別（&：かつ）
if (a > 0 || a < 2) { "a > 0 or a < 2" } else { "no" } # a > 0またはa < 2 か否かを判別（||：または）
if (a > 0) { "positive" } # 0より大の時だけ動作、否の場合は何も起こらない（空白）
if (a > 0 || a < 2) {
"a > 0 or a < 2"
} else {
"no"
} # 条件式が長い場合など、改行してもOK
# if ... 、複数の if, else を重ねることもできます（複合的な条件分岐）。
# 条件判断と返す値の部分が、複数の場合。上で定義した、
f <- c(2:6) # を例に
if (f > 3) { "over 3" } else { "under 3" } # if, elseでは、fの最初の一つ目の値で条件判断する
# fの各値について条件判断したい場合はifelseを使用する
f <- ifelse(f > 3, f*2, 0) # 3より大きければ２倍、小さければ0
# 注：ifelseは、条件式と返す値の個数を対応させる必要あり（個数が同じベクトル）。
# 個数が異なる場合、条件式の形に合う部分しか値を返さない（警告もしてくれない）
##### 繰り返し作業の自動化＝プログラミングの初歩 #############
sapply(1:5, function(s) sample(1:5, 10, replace=T))
# 関数を定義しつつ、繰り返し計算の結果を行列に結合するなど、荒っぽい処理でも対応しやすい
# ただし帰ってくる値が、ちゃんと目的のものに合っているかどうか、
# 結果が明確なよりわかりやすい例でチェックすることを習慣づけるべきです
# もうひとつの繰り返し作業の例（for ループという）
# 他の言語でもよく使われるが、ちょっと処理速度は遅め
x <- 0 # xの元の値は0
for (i in 1:5) { x <- x + 1 } # 繰り返し範囲の始点と終点に {} を付ける
# iを1~5まで繰り返し。左辺のxは新しく定義するx、右辺のは直前に定義したx
x # 繰り返し作業後のxを入力すると値が変わっているのが確認できる
# 7行おきに合計数を集計
xx <- matrix(1:56,28,2) # 集計をする元データ（28行、2列）
xxx <- matrix(0,4,2) # 集計結果を納める入れ物（列毎に集計）
for (i in 1:4) {
xxx[i,1] <- sum(xx[((7*(i-1)+1):(7*i)),1]) # [ ] 内の式は、なぜか掛算でも( ) にしないと計算が狂う
xxx[i,2] <- sum(xx[((7*(i-1)+1):(7*i)),2])
}
xxx
# 複数のファイルを連続自動読み込み
# あるディレクトリ（フォルダ）にあるcsvファイルを全部読み込む場合
data.lists <- list() # データを格納する入れ物
files <- list.files() # ファイルリスト（csv以外のファイルが混じってると次でエラーになる）
for (i in 1:(length(files))) { # filesの個数分、読み込む
data.lists[[i]] <- read.csv(file=files[i]) }
data.lists[[1]] # 1番目のファイル内容を呼び出す場合
##### RやRの機能を論文に引用する場合 #####################
citation() # 使用しているRを引用する場合のreference情報
citation("パッケージ名") # 特定のパッケージのreference情報
citation("stats") # 関数 glm を引用する場合
citation("lme4") # 関数 lmer を引用する場合
0 件のコメント:
コメントを投稿