最初の方でやること
data <- read.delim("clipboard", header=T) # Excel からコピー&ペースト（for Win）。列名を入れたい場合はheader=T。
data <- read.delim(pipe("pbpaste"), header=T) # Excel からコピー&ペースト（for Mac）
head(data) # 中身確認
se <- function(x){sd(x)/sqrt(length(x));} # Standard Error 計算式
d <- subset(d, treatment=="Negative") #データフレームからある列に対して条件を満たした行のみ抽出
d$treatment <- factor(d$treatment, levels=c("Control", "1ng", "10ng") #統計やグラフで処理区名の順を任意に出力するためにあらかじめ並び替えておく
作成したテーブルを保存（ディレクトリ先のフォルダに保存）
write.table(x, "Relative abundance of each compound.csv", sep=",") # 変数 x に格納した表をCSV形式（sep=","でコンマ区切り）で書き出し保存
write.csv(x, "Data of termite chemicals.csv", row.names=F) #変数xに格納したデータフレーム等をCSVファイルとして書き出し保存。row.namesをFALSEにしておくと１列目が行番号にならずに済む。
その他に役立つコマンド
for(i in 1:20){
x <- 100 + i
variable <- paste("t", i, sep="")
assign(variable, x)
}
#for文の中で出力された結果 x（数、データフレーム等。上の例では100に i を足した結果）をそれぞれ別の変数（ここではt1 ~ t20と番号付けした）に 格納したい場合
下準備
library(dplyr)
library(reshape2)
library(ggplot2)
dataMean <- tapply(d[, 9], d[, c(2, 4)], mean) #平均を計算
dataSe <- tapply(d[, 9], d[, c(2, 4)], se) #標準誤差を計算
dmean = reshape2::melt(dataMean, value.name="mean") # ggplot2用に平均について行列並べ替え
dse = reshape2::melt(dataSe, value.name="se") # ggplot2用にSEについて行列並べ替え
g <- merge(dmean, dse) #平均とSEの一覧表を一つの行列にまとめる
g
str(g) # gにどのようなデータが入っているか確認
棒グラフ（縦向き）
p <- ggplot(g, aes(x=Sex, y = mean)) #ggplot2の設定
p <- p + geom_bar(stat="identity", width=0.6) # 棒グラフの設定
p <- p + theme_classic() #グラフのデザインの設定（classic：白黒、ガイド線無）
p <- p + scale_y_continuous(breaks=seq(0, 15, by=2)) #縦軸の目盛の刻み方の設定（この場合、範囲は0~15で2刻みに目盛）
p <- p + theme(axis.title = element_text(size=8),
axis.text.x = element_text(angle=45, hjust=1, size=8),
axis.text.y = element_text(size=8)) #縦軸・横軸タイトルの文字サイズと横軸ラベル45°回転
p <- p + geom_errorbar(aes(ymin = mean - se, ymax = mean + se, width = 0.3)) #エラーバー描画
p <- p + labs(x="", y="Abundance (Mean ± SE) [μg]") #軸ラベル挿入
p <- p + facet_wrap(~Soldier) #グラフをカテゴリごとに分割（カテゴリ毎に縦軸の範囲変える時は facet_wrap(~ Soldier, scales="free_y")）
dev.new() #新しいウィンドウを表示 括弧の中に(width=5, height=5)と書けばウィンドウの寸法指定可能（単位：inch）
plot(p) #グラフ描画
左右に伸びる棒グラフ（使用用途：Two-choice testの結果など）
library(gridExtra)
# Mean and SEM of prop. of workers on left papers
dataLmean <- tapply(d[, 9], d[, 3], mean)
dataLse <- tapply(d[, 9], d[, 3], se)
dLmean <- melt(dataLmean, value.name="mean")
dLse <- melt(dataLse, value.name="se")
becL <- merge(dLmean, dLse)
Paper <- c("Left", "Left", "Left")
becL <- cbind(Paper, becL)
# Mean and SEM of prop. of workers on right papers
dataRmean <- tapply(d[, 10], d[, 3], mean)
dataRse <- tapply(d[, 10], d[, 3], se)
dRmean <- melt(dataRmean, value.name="mean")
dRse <- melt(dataRse, value.name="se")
becR <- merge(dRmean, dRse)
Paper <- c("Right", "Right", "Right")
becR <- cbind(Paper, becR)
#処理区名の部分を生成
gbec_tr <- ggplot(data=becR, aes(x=1,y=Var1)) + geom_text(aes(label=Var1),size=3) + ylab(NULL) + ggtitle("") +
scale_x_continuous(expand=c(0,0),limits=c(1,1)) +
scale_y_discrete(limits=c("XXX vs YYY", "XXX vs ZZZ", "Hexane vs Hexane")) + theme_classic()+
theme(axis.title=element_blank(), panel.grid=element_blank(), plot.title = element_text(size = 9),
axis.text.y=element_blank(), axis.ticks.y=element_blank(),
axis.line.x = element_blank(), axis.line.y = element_blank(), panel.background=element_blank(),
axis.text.x=element_text(color=NA), axis.ticks.x=element_line(color=NA), plot.margin = unit(c(2, 0, 3, 0), "mm"))
# Left 側を生成
gbec_left <- ggplot(data=becL, aes(x=Var1, y=mean)) + geom_bar(stat="identity", width=0.8) + ggtitle("") +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se, width = 0.3)) +
scale_x_discrete(limits=c("XXX vs YYY", "XXX vs ZZZ", "Hexane vs Hexane")) +
scale_y_reverse(expand=c(0, 0), limits=c(0.64, 0), breaks=c(0.6, 0.4, 0.2)) + coord_flip() + theme_classic()+
theme(axis.title.x = element_blank(), axis.title.y = element_blank(),
plot.title = element_text(size = 9), axis.line.y = element_blank(), axis.text.y = element_blank(),
axis.ticks.y = element_blank(), plot.margin = unit(c(2, -0.9, 3, 0), "mm"))
# Right 側を生成
gbec_right <- ggplot(data=becR, aes(x=Var1, y=mean)) + geom_bar(stat="identity", width=0.8) + ggtitle("") +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se, width = 0.3)) +
scale_x_discrete(limits=c("XXX vs YYY", "XXX vs ZZZ", "Hexane vs Hexane")) +
scale_y_continuous(expand=c(0, 0), limits=c(0, 0.64)) + coord_flip() +theme_classic()+
theme(axis.title.x = element_blank(), axis.title.y = element_blank(), plot.title = element_text(size = 9),
axis.line.y = element_blank(), axis.text.y = element_blank(),
axis.ticks.y = element_blank(), plot.margin = unit(c(2, 1, 3, -0.9), "mm"))
# 左から処理区名部分・Left側・Right側の順に組み合わせる
dev.new(width=4, height=5)
grid.arrange(gbec_tr, gbec_left, gbec_right, ncol=3, widths=c(3/10, 3.5/10, 3.5/10))
折れ線グラフ
p <- ggplot(g1, aes(x=Var1, y = mean, group=1)) #折れ線が1本の時はgroup=1、グループごとに複数本描く場合はgroup=（グループ名） ＊折れ線が1本の時にgroupを指定しないと線が描画されない。
p <- p + geom_point() + geom_line() #折れ線グラフの各点とそれをつなぐ線
p <- p + ylim(0, 2.5) #Y軸の上限・下限を指定
p <- p + theme_classic() #グラフのテーマ指定（この場合は最もシンプルなclassic）
p <- p + scale_x_discrete(limits = c("Negative", "(1ng)", "(10ng)", "(100ng)"), labels=c("(1ng)"="1 ng", "(10ng)"="10 ng", "(100ng)"="100 ng")) #limitsで横軸のカテゴリの順序変更、labelsで名称変更
p <- p + geom_errorbar(aes(ymin = mean - se, ymax = mean + se, width = 0.5)) #エラーバー描画
p <- p + labs(x="", y="Antennation (Mean ± SE)", title="Sample") + theme(axis.title = element_text(size=9), axis.text.x = element_text(angle=45, hjust=1, size=9), axis.text.y = element_text(size=9)) #labsで軸ラベルやタイトル挿入、themeで書式設定
dev.new() #新しいグラフ描画ウィンドウを出す
plot(p) #グラフ描画ウィンドウにグラフを描画
回帰直線付き
p <- ggplot(d, aes(x=volume, y=weight))
p <- p + geom_smooth(method="lm", colour="black", size=0.5) + geom_point() + theme_classic() +
labs(x="Total volume [mL]", y="Total weight [g]") +
theme(axis.title = element_text(size=8), axis.text = element_text(size=7))
dev.new(width=4, height=3)
plot(p)
# geom_smooth(method="lm", colour="black", size=0.5)で回帰直線を引ける。methodで回帰式の種類を変更可。lm、glm、nls、rlm、loess、gamを指定可。glmの場合、family=""で分布を指定可能。fomulaで式を指定できる。se=FALSEで信頼区間を非表示に。levelで信頼区間を変更可（デフォルトは95%信頼区間。level=0.95）。
多項回帰の場合はmethod="lm"で、fomula = y ~ x + I(x^2) + I(x^3) +…と、2次以上の項はI( )で囲う必要あり
PCA (Principle component analysis 主成分分析)
rpca <- prcomp(dp) # dpには行に各サンプル、列に変数の並んだデータフレーム
summary(rpca) #第1、第2、…第n主成分の重要度と累積重要度を表示
rpca$rotation #各主成分に対するそれぞれの変数の寄与率表示
# プロット作成
library(ggfortify)
autoplot(prcomp(dp[, 3:15]), data=dp, shape="Caste", colour="Colony", loadings=T, loadings.label=T, loadings.label.size=3) + scale_shape_manual(values=c(16,7,8,15,17,18)) + theme_bw()
# shapeの種類はRのデフォルトのグラフィックパラメーターに準拠
線形判別分析 (LDA: Linear discriminant analysis)
library(MASS)
library(ggplot2)
z <- lda(Sample~., data=d)
z
table(d[,2], predict(z)$class) #判別結果
data.lda <- data.frame(No=c(1:7), coef(z)) #Extract scaling for each predictor (top 7 predictors)
data.lda$length <- with(data.lda, sqrt(LD1^2 + LD2^2))
dd <- cbind(d[, 1:4], predict(z)$x)
p <- ggplot(dd, aes(x=LD1, y=LD2, colour=Sample, shape=Colony)) #固有ベクトル入りのプロットを描画
p <- p + geom_point(size=3) + theme_bw() + scale_shape_manual(values=c(15:19, 7)) +
geom_text(data=data.lda, aes(x=LD1*0.009, y=LD2*0.009, label=No, shape=NULL, linetype=NULL), size = 3, vjust=0.5, hjust=0, color="red") +
geom_segment(data=data.lda, aes(x=0, y=0, xend=LD1*0.0085, yend=LD2*0.0085, shape=NULL, linetype=NULL), arrow=arrow(length=unit(0.2,"cm")), color="red")
dev.new(width=6.4, height=5)
plot(p)
# geom_text()内は固有ベクトルの成分名の文字のサイズ・色・配置を設定。geom_segment()内で固有ベクトルの色や（相対）長さを設定
マススペクトル描画用
library(ggplot2)
p <- ggplot(d, aes(x=m.z, y=relative_intensity)) # dにはm/z値と相対強度 (%) の2列を含むdata frameを格納
p <- p + geom_bar(stat="identity", width=0.6) + theme_classic()
p <- p + scale_y_continuous(breaks=seq(0, 100, by=25), expand=c(0,0)) # expandで横軸とbarの間の隙間調整
p <- p + xlim(10, 290) + scale_x_continuous(breaks=seq(10, 290, by=10))
p <- p + labs(x="m/z", y="Relative intensity [%]", title="")
plot(p)
t test
t.test(y ~ x, data = d, var.equal=T) # データフレーム d において従属変数 y ついて、独立変数 x に含まれる二群間で比較
t.test(a, b, var.equal=T) #ベクトルaのグループとベクトルbのグループ間で比較
# var.equal が T …対応のあるt検定、F…対応のないt検定（Welchのt検定）
Wilcoxonの順位和検定（＝Mann-Whitney's U test）
wilcox.test(a, b) #ベクトルaのグループとベクトルbの二群間で比較
pairwise.wilcox.test(x, g, p.adjust.method = bonferroni) #xの値についてgに格納されたグループ間で多重比較
正規性の検定
shapiro.test(data$area) #Shapiro-wilk検定
tapply(data$area, data$caste, shapiro.test) #カテゴリ毎にShapiro-wilk検定
等分散性の検定
var.test(a, b) #ベクトルaのグループとベクトルbの二群間で等分散検定（F検定）
leveneTest(Area ~ Treatment, data=d) #正規母集団以外でも適用可能な等分散検定。多群比較も可（ルビーン検定）
二項検定 Binomial test
binom.test(c(n1, n2)) # n1, n2 はカウント数。(n1 + n2)に対するn1の割合とn2の割合を比較
binom.test(c(n1, n2), p = 1/2) #コインを多数回投げた時に表と裏の出る確率が等しいか検定する場合（n1: 表が出た回数、n2: 裏が出た回数）
分散分析 Analysis of variance (ANOVA)
r <- lm(y ~ x, data=d) #一元配置分散分析をする場合
r <- lm(y ~ x1 + x2, data=d) #二元配置分散分析をする場合（交互作用見ないなら"+"、見るなら" * "）
summary(aov(r)) #結果を出力
TukeyHSD(aov(r)) #続けてTukey's HSD testを行う場合
# multcompライブラリを使うと多重比較した後、以下のように自動でアルファベットをつけることもできる
library(multcomp)
r <- lm(y ~ x, data=d)
res <- aov(r)
res2 <- glht(res, linfct=mcp(x="TukeyHSD"))
cld(res2, level=0.05, decreasing=T) # 有意水準0.05、降順でアルファベットをつける場合
回帰分析 Regression analysis
r <- lm(y ~ x, data=d) #データフレームdにある独立変数xに対する従属変数yの単回帰
summary(r) #結果出力
順序ロジスティック回帰
library(ordinal)
r <- clm(as.factor(Area) ~ Colony*Treatment, data=d) #モデルを指定（従属変数はFactor型である必要がある）
anova(r, type="II")
# 多重比較する場合
library(emmeans)
library(multcomp)
result1 <- emmeans(r, ~ Treatment + Colony) #左のようにかけばTreatmentとColonyの全組み合わせで比較可能
pairs(result1, adjust="tukey") #Tukey's HSD testで多重比較
cld(result1, Letters=letters) #有意に異なる組み合わせで異なるアルファベットを自動でつけてくれる
一般化線形混合モデル Generalized Linear Mixed Model (GLMM)
library(lme4)
library(multcomp)
r <- glmer(Abundance ~ Treatment + (1|Colony), family="poisson", data=d) # ポアソン分布の場合
r <- glmer(cbind(sample, total-sample) ~ Treatment + (1|Colony), family="binomial", data=d) #二項分布の場合
# Convergence warningが出た場合 "glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 10000))" の項を入れる
r <- lmer(Abundance ~ Treatment + (1|Colony), data=d) # 正規分布の場合＝線形混合モデル Linear Mixed Model (LMM)
rNULL <- lmer(Abundance ~ 1 + (1|Colony), data=d) #対立モデル
anova(r, rNULL, test="Chisq") # 分散分析
rT <- glht(r, linfct=mcp(Treatment="Tukey")) # 続けて多重比較（Tukey HSD test）を行う場合
rT <- glht(r, linfct=mcp(Treatment="Dunnett")) # 続けて多重比較（Dunnett's test の場合）を行う場合
summary(rT) # 結果を出力
cld(rT, level=0.05, decreasing=T) #降順（decreasing=Tの場合）でアルファベットをつける
多変量分散分析 Multivariate Analysis of Variance (MANOVA)
result1 <- manova(cbind(dp[, 3], dp[, 4], dp[, 5], dp[, 6],dp[, 7], dp[, 8], dp[, 9], dp[, 10], dp[, 11], dp[, 12], dp[, 13], dp[, 14], dp[, 15]) ~ Caste*Sex, data=dp) #従属変数が多数あるのでcbindでまとめる
result1 #要約統計量を表示
summary(result1, test="Pillai", tol=0) #Pillaiのトレース検定
summary(result1, test="Wilks", tol=0) #Wilk's lambda test
summary(result1, test="Hotelling-Lawley", tol=0) #Hotelling-Lowleyのトレース検定
summary(result1, test="Roy", tol=0) #Royの最大根検定
summary.aov(result1) #各変数ごとにANOVA。表示されるResponse No.が各従属変数に該当
Permutational Multivariate Analysis of Variance (PERMANOVA)
library(vegan)
result1 <- adonis(cbind(dp[, 3], dp[, 4], dp[, 5], dp[, 6],dp[, 7], dp[, 8], dp[, 9], dp[, 10], dp[, 11], dp[, 12], dp[, 13], dp[, 14], dp[, 15]) ~ Caste, permutations=1000, data=dp, method="bray") #従属変数が多数あるのでcbindでまとめる
result1 #統計結果を表示
Fisher's exact test
# 2群間比較の場合
mx <- matrix(c(2, 10, 15, 3), nrow=2, byrow=T) #分割表を用意（Matrix型でもデータフレーム型でも可）
fisher.test(mx)
# 3群以上の比較の場合
pairwise.fisher.test(a$Sample, a$Total, p.adjust.method="bonferroni") #a はデータフレームも可（この場合aの各行が比較した処理群）。Bonferroni補正で比較する場合は"bonferroni"と書く