必要なパッケージ
require(here)
require(plotly)Rとは
Rは統計解析のためのフリーソフトです。もう少し正確に言うと、Rはコンピュータ言語の名前であり、パソコンにソフトウェアとしてインストールされたRは、R言語を使用するための「環境」です。Rは多くの機能を持ち、統計解析からデータの前処理、データの概観、論文用グラフの作成まで幅広く利用することができます。また、「パッケージ」として配布されている拡張プログラムをインストールすることで、様々な解析を簡単に行うことができます。このため、Rを使いこなすスキルは、農学や生命科学の研究者にとって非常に有用なものとなっています。
Rを使った簡単な計算
Rを使った解析では、基本的にコマンドを逐次入力しながら対話的に解析を進めていきます。ただし、実際に解析を行う際には、事前に一連のコマンドをRスクリプトとして用意しておくと便利です。Rスクリプトを実行することで、部分的な修正ができたり、解析の履歴を確認できたりするため、これを覚えておくとよいでしょう。
まずは、コマンド入力で簡単な計算を行いながら、Rに慣れることから始めます。
Rの最も簡単な使い方は、簡単な算術式を入力して答えを出すことです。例えば
3 + 5 * 3## [1] 18得られた結果をもとに次の計算を行いたい場合は、次のように何らかの変数に結果を代入します。
x <- 1 + 2
x## [1] 3代入された値は、作成した変数の名前を通して別の計算に使用することができます。
x + 5 * x## [1] 18関数を用いて様々な計算を行うことができます。
abs(x) # 絶対値## [1] 3sin(x) # 正弦## [1] 0.14112log(x) # 底eの対数## [1] 1.098612exp(x) # 指数## [1] 20.08554factorial(x) # 階乗## [1] 6もう少し複雑な計算をしてみましょう。平均が \(mu\) で分散が \(sigma^2\)の正規分布の確率密度関数は、次のようになります。 \[f(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\] これをRを使って計算してみましょう。
mu <- 3
s2 <- 2
x <- 5
1 / sqrt(2 * pi * s2) * exp(- (x - mu)^2 / (2 * s2)) ## [1] 0.1037769計算結果の確認のために、dnorm関数で正規分布の確率密度を計算しても、同じ結果が得られます。
dnorm(x, mu, sqrt(s2))## [1] 0.1037769
図1. 平均3、分散2の正規分布
確認クイズ
次の分布はポアソン分布と呼ばれる分布で、起こりにくい事象が生じる回数の確率などに用いられます。起こりにくい事象が生じる回数の平均が\(\lambda\)のときに、その事象が\(x\)回生じる確率は、 \[ P(x = k) = \frac{\exp(-\lambda) \lambda ^k}{k!} \] と表せます。この式を用いて、\(\lambda = 0.3\)のとき、その事象が3回生じる（\(x = 5\)）確率を計算してみましょう。
lmbd <- 0.3
x <- 4
# 自分でコードを書いてみましょう。計算結果の確認は、dpois関数を用いて行うことができます。
dpois(x, lmbd)## [1] 0.0002500261簡単な統計量の計算
まず、8個の要素からなるベクトルを作ってみましょう。これは、ブルーベリーの糖度（可溶性固形物含量、soluble solid content: SSC）のデータです（データの出典は後述します）。
sug <- c(10.0, 10.5, 13.0, 11.9, 10.4, 8.6, 10.0, 10.2)
sug## [1] 10.0 10.5 13.0 11.9 10.4 8.6 10.0 10.2同じサンプルの酸度（滴定酸度、titratable acidity: TA）も入力して、糖度と酸度の比をとってみましょう。
aci <- c(0.797, 0.192, 0.219, 0.839, 0.511, 0.449, 0.162, 0.158)
rat <- sug / aci
rat## [1] 12.54705 54.68750 59.36073 14.18355 20.35225 19.15367 61.72840 64.55696期待値（平均） \[ m_x = \frac{\sum_{i=1}^n x_i}{n} \] は、mean関数を用いて計算できます。
mean(rat)## [1] 38.32126分散の推定値 \[ v_x = \frac{\sum_{i=1}^n (x_i - m_x)^2}{n - 1} \] は、var関数を用いて計算できます。
var(rat)## [1] 554.84852変数間の共分散の推定値 \[ v_{x,y} = \frac{\sum_{i=1}^n (x_i - m_x)(y_i - m_y)}{n - 1} \] は cov関数で計算できます（ここで、\(m_y\)は変数\(y\)の期待値）。
ここでは、このブルーベリーの官能評価をした結果を新たな変数として入力し、糖度と酸度の比との共分散を計算してみましょう。
lik <- c(19.7, 25.0, 28.2, 16.1, 23.5, 18.9, 21.9, 27.3)
cov(rat, lik)## [1] 78.63015相関係数 \[ r_{x,y}=\frac{\sum_{i=1}^n (x_i - m_x)(y_i - m_y)}{\sqrt{\sum_{i=1}^n (x_i - m_x)^2}\sqrt{\sum_{i=1}^n (y_i - m_y)^2}} \] は、cor関数を使って計算できます。
cor(rat, lik)## [1] 0.7904019官能評価の結果と、糖度と酸度の比には、比較的高い相関があることがわかります。
なお、cor.test関数を使うと、検定して統計的に有意であるか確認することもできます。
cor.test(rat, lik)##
## Pearson's product-moment correlation
##
## data: rat and lik
## t = 3.1605, df = 6, p-value = 0.01955
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1935081 0.9602434
## sample estimates:
## cor
## 0.7904019\(p\)値は、5%水準で有意です。
なお、今回の講義では、時間的な問題から、統計的に非常に重要な仮説検定については、説明をしません。この講義では、データマイニングなどとよばれる手法（収集された情報のなかから傾向や関連性を見出す分析手法）について説明します。しかし、基礎的な統計学は、データマイニングを行うためにも重要な知識です。是非、別の機会にしっかりと勉強してください。
外部データの取り込みと解析
Rを自分の研究に使う場合、表計算ソフトで整理したデータを読み込んで分析することがほとんどです。ここでは、データとして、[Gilbertら（2015, PLoS ONE 10: e0138494）] (https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0138494) を用います。この研究は、複数のブルーベリー品種について多環境試験（複数試験地、複数年次に渡る栽培試験）を行い、嗜好性などの官能評価の結果と、酸度、糖含量、多数の揮発性成分などの生化学的特性との関係について解析を行ったものです。 データは、同論文のSupporting Informationに含まれています。
ここでは、上述したデータから、すべての環境で栽培されていた５品種について揮発性成分以外のデータを抜き出し、一部細工をして（後述します）csvファイルとして保存したデータを準備しました。
bb <- read.csv(here("data", "blueberry_5var_miss.csv"), stringsAsFactors = T) # csvファイルの読み込み
# here関数に、データの格納されているフォルダとファイル名を入れると、アクセスできる
# stringsAsFactors = Tは、文字列をfactorとよばれる形式で読み込むように指定している取り込んだデータのサイズやデータの一部を確認するには、次のようにします。
dim(bb)## [1] 133 19head(bb)## gID Year Location Harvest Genotype Panelists Overall.Liking
## 1 2012 C1 Emerald 2012 C 1 Emerald 95 23.5
## 2 2012 C1 Endura 2012 C 1 Endura 95 17.5
## 3 2012 C1 Farthing 2012 C 1 Farthing 95 24.2
## 4 2012 C1 Meadowlark 2012 C 1 Meadowlark 95 22.3
## 5 2012 C1 Primadonna 2012 C 1 Primadonna 95 24.0
## 6 2012 C2 Emerald 2012 C 2 Emerald 90 20.0
## Texture Sweetness Sourness Flavor SSC TA SSC.TA pH Fructose Glucose
## 1 24.9 22.4 10.7 25.0 10.4 0.5110 20.4 3.43 48.74 49.39
## 2 25.1 19.2 22.6 27.8 10.5 0.6199 16.9 3.36 47.79 46.49
## 3 23.9 23.6 11.8 26.5 10.6 0.3064 34.6 3.55 51.31 51.66
## 4 28.7 21.5 11.2 24.4 9.2 0.2631 34.8 3.78 42.57 45.11
## 5 23.6 25.5 10.4 25.8 11.3 0.2346 48.2 3.95 56.83 56.35
## 6 23.0 20.2 13.6 22.3 10.6 0.3992 26.6 3.54 48.12 49.15
## Sucrose Total.Sugar
## 1 0.13 98.26
## 2 0.19 94.47
## 3 0.18 103.15
## 4 0.12 87.79
## 5 0.21 113.39
## 6 0.22 97.48データの分析
多数の変数について、その分布や変数間の関係を調べるために解析を行いたいことがよくあります。しかし、計測データには、実験的な理由による欠損項目が含まれている場合があります。ここでは、先にcsvファイルから読み込んだデータを解析することで、その対策を考えます。
先程と同じ方法で、糖度と酸度の比をとり、その平均値を計算してみます。
rat <- bb$SSC / bb$TA
mean(rat)## [1] NA平均値を算出することができず、“NA”という値しか得られません。なぜでしょうか？
これは、糖度や酸度に欠損値（RではNAと表現）が含まれているためです。
rat## [1] 20.352250 16.938216 34.595300 34.967693 48.167093 26.553106
## [7] 18.141429 32.095097 37.509377 57.291667 33.252721 32.840722
## [13] 48.021343 61.652281 76.000000 12.554928 14.899957 NA
## [19] 19.153675 20.329881 15.543329 16.734439 13.003599 26.851098
## [25] 23.399476 20.591341 15.189873 14.395140 45.131634 34.620174
## [31] 9.544116 10.517892 43.936731 11.800259 19.379149 11.518822
## [37] 19.218025 49.335863 18.947793 29.907429 13.737494 24.688562
## [43] NA 27.168879 26.669990 17.042688 31.257691 20.285935
## [49] 54.659032 72.134962 28.267635 30.069391 29.049055 49.648947
## [55] 89.852827 32.743067 52.754435 47.722343 44.483210 58.350101
## [61] 15.494031 13.130445 14.082704 15.848171 28.419183 14.401935
## [67] 11.518546 14.454026 16.910739 25.261721 34.261838 22.519509
## [73] 21.693625 22.871665 15.274529 18.449497 13.544669 25.668166
## [79] 28.637177 28.274968 15.384615 16.812609 30.754289 NA
## [85] 36.823105 14.751175 40.702314 53.878034 36.654804 38.759690
## [91] 51.756007 64.458900 64.638783 NA 59.976932 57.917109
## [97] 82.797428 145.234493 25.503356 22.176591 41.325696 77.447336
## [103] 68.027211 22.749559 28.289298 31.127314 19.306541 62.569214
## [109] 19.004671 20.912893 56.898288 34.727144 54.642381 30.855916
## [115] 30.380466 30.835735 52.322738 NA 23.956723 38.167939
## [121] 21.938192 55.045872 59.306569 14.176793 30.326005 32.167475
## [127] 41.133101 44.226044 27.549541 47.148630 56.338028 NA
## [133] 70.101857対象とする変数に欠測値が含まれて計算結果が得られない場合は、計算時にオプションとしてna.rm = Tを指定してください。
mean(rat, na.rm = T)## [1] 34.23757計算されるようになりました。
次に、官能試験の結果と糖度と酸度の比の相関係数を計算してみましょう。
lik <- bb$Overall.Liking
cor(rat, lik)## [1] NA結果はNAとなります。これは、前回と同様にデータの欠損が原因です。
欠損値を処理するオプションを指定し、再度計算を試みます。相関計数を計算するcorや、共分散を計算するcovなど、2変数以上の関係を表す統計量を計算するときは、欠測値に対する値をna.rm = Tではなくuse =というオプションで指定することに注意すます。これは、変数A, B, Cがあった場合、それぞれの変数について異なるサンプルが欠測している場合があるため、欠測の扱いが、欠測を無視する（na.rm = T）という単純な指定だけでは済まなくなるためです。これを確認してみましょう。
cor(bb$Overall.Liking, rat, use = "pair")## [1] 0.3975358今回は、相関関係の算出に成功します。
次に、官能評価の結果、糖度、酸度、の3変数の間にどのような相関があるか確認してみます。
df <- data.frame(lik = bb$Overall.Liking, sug = bb$SSC, aci = bb$TA)
cor(df)## lik sug aci
## lik 1 NA NA
## sug NA 1 NA
## aci NA NA 1欠測のために、相関係数は、やはり計算されません。
cor(df, use = "pair")## lik sug aci
## lik 1.0000000 0.40315843 -0.47418923
## sug 0.4031584 1.00000000 -0.05795763
## aci -0.4741892 -0.05795763 1.00000000今度は、計算されました。官能評価と、糖度、酸度には、それぞれ正、負の相関があることがわかります。また、酸度と糖度には大きな相関はみられません。
確認クイズ
欠測の扱いかたを、別のオプションを使って、変更してみましょう。
use = “pair”は、use = “pairwise.complete.obs”の略で、2変数間で欠測がないサンプルをもとに相関係数が計算されます。一方、use = “comp”は、use = “complete.obs”の略で、対象としている変数（ここでは3つの変数）のすべてで欠測がないサンプルだけで相関係数が計算されます。
# 自分でコードを書いてみましょう。計算される相関係数を前のものと比較してみてください。
データの再読み込みと変数の設定
実は、糖度と酸度に含まれていた欠測値は、オリジナルのデータでは、計測されているものでした。そこで、改めて、欠測のないデータを読み込みます。
bb <- read.csv(here("data", "blueberry_5var.csv"), stringsAsFactors = T)データの内容を確認します。
colnames(bb) # 列名の表示## [1] "gID" "Year" "Location" "Harvest"
## [5] "Genotype" "Panelists" "Overall.Liking" "Texture"
## [9] "Sweetness" "Sourness" "Flavor" "SSC"
## [13] "TA" "SSC.TA" "pH" "Fructose"
## [17] "Glucose" "Sucrose" "Total.Sugar"変数を要約して確認します。
summary(bb) # summary関数で、値の分布などを確認できる## gID Year Location Harvest
## 2012 C1 Emerald : 1 Min. :2012 C:44 Min. :1.000
## 2012 C1 Endura : 1 1st Qu.:2012 H:44 1st Qu.:1.000
## 2012 C1 Farthing : 1 Median :2013 W:45 Median :2.000
## 2012 C1 Meadowlark: 1 Mean :2013 Mean :1.992
## 2012 C1 Primadonna: 1 3rd Qu.:2014 3rd Qu.:3.000
## 2012 C2 Emerald : 1 Max. :2014 Max. :3.000
## (Other) :127
## Genotype Panelists Overall.Liking Texture
## Emerald :27 Min. : 72.00 Min. :12.40 Min. :12.9
## Endura :27 1st Qu.: 87.00 1st Qu.:18.90 1st Qu.:22.5
## Farthing :27 Median : 91.00 Median :21.90 Median :24.7
## Meadowlark:27 Mean : 90.83 Mean :21.81 Mean :24.3
## Primadonna:25 3rd Qu.: 95.00 3rd Qu.:24.90 3rd Qu.:26.8
## Max. :108.00 Max. :30.00 Max. :31.5
##
## Sweetness Sourness Flavor SSC
## Min. :16.30 Min. : 5.90 Min. :17.00 Min. : 8.10
## 1st Qu.:20.20 1st Qu.:10.90 1st Qu.:24.00 1st Qu.:10.00
## Median :22.60 Median :14.00 Median :25.60 Median :10.70
## Mean :22.52 Mean :15.15 Mean :25.54 Mean :10.75
## 3rd Qu.:25.00 3rd Qu.:19.00 3rd Qu.:27.20 3rd Qu.:11.40
## Max. :30.70 Max. :28.10 Max. :31.40 Max. :13.90
##
## TA SSC.TA pH Fructose
## Min. :0.0661 Min. : 9.5 Min. :2.820 Min. :28.40
## 1st Qu.:0.2293 1st Qu.: 19.2 1st Qu.:3.230 1st Qu.:43.40
## Median :0.3779 Median : 30.1 Median :3.490 Median :47.72
## Mean :0.4155 Mean : 34.6 Mean :3.489 Mean :47.78
## 3rd Qu.:0.5402 3rd Qu.: 47.0 3rd Qu.:3.710 3rd Qu.:51.93
## Max. :0.9849 Max. :145.2 Max. :4.410 Max. :69.86
##
## Glucose Sucrose Total.Sugar
## Min. :31.48 Min. :0.040 Min. : 60.66
## 1st Qu.:44.01 1st Qu.:0.240 1st Qu.: 89.90
## Median :48.36 Median :0.840 Median : 97.07
## Mean :49.02 Mean :1.432 Mean : 98.25
## 3rd Qu.:52.63 3rd Qu.:2.060 3rd Qu.:107.05
## Max. :76.01 Max. :5.440 Max. :140.36
## 年次（Year）、収穫時期（Harvest）、パネリスト（Oanelists）のIDが数値データとなっているので、factor（因子）に変換します。
bb$Year <- as.factor(bb$Year)
bb$Harvest <- as.factor(bb$Harvest)
bb$Panelists <- as.factor(bb$Panelists)
summary(bb[,1:8])## gID Year Location Harvest Genotype
## 2012 C1 Emerald : 1 2012:45 C:44 1:45 Emerald :27
## 2012 C1 Endura : 1 2013:44 H:44 2:44 Endura :27
## 2012 C1 Farthing : 1 2014:44 W:45 3:44 Farthing :27
## 2012 C1 Meadowlark: 1 Meadowlark:27
## 2012 C1 Primadonna: 1 Primadonna:25
## 2012 C2 Emerald : 1
## (Other) :127
## Panelists Overall.Liking Texture
## 90 :25 Min. :12.40 Min. :12.9
## 95 :24 1st Qu.:18.90 1st Qu.:22.5
## 93 :15 Median :21.90 Median :24.7
## 85 :10 Mean :21.81 Mean :24.3
## 101 :10 3rd Qu.:24.90 3rd Qu.:26.8
## 72 : 5 Max. :30.00 Max. :31.5
## (Other):443年間4試験地で3回収穫されて計測されたデータであることが分かります。
データの可視化
実際に統計解析を行う前に、データを様々な角度から見てみることは非常に重要です。例えば、前述した平均や分散といった統計は要約のための統計であり、平均や分散が似ている変数であっても、観測値の分布が大きく異なる場合があります。そのため、データをよく見て、データの特徴を把握することが重要です。また、分析結果を論文として発表するための準備をする際にも、データの可視化が必要となります。ここでは、様々なデータの可視化手法を学びます。
まず、可視化手法の説明の前に、データ中の変数を直接呼び出せるようにします。
attach(bb)attachコマンドにより、bb$…という記述をしなくても変数を指定できるようになります。このコマンドを使わない場合は、bb$変数名というかたちで変数を指定しなければなりません。
まずは、糖度のヒストグラムを描いてみます。
hist(SSC)
次に、箱ひげ図を描いてみます。
boxplot(SSC)
次に、段階評価のデータの図示について説明するために、官能試験の結果であるOverall.Likingについて「無理やり」10段階評価に変換してみます。
score <- (Overall.Liking - min(Overall.Liking)) /
(max(Overall.Liking) - min(Overall.Liking)) # 0から1の値をとるスコアにする
score <- floor(score / (1 + 1e-10) * 10)
# 1未満になるように1よりごくわずか大きくなる値で除算して10倍し、小数点以下を切り捨てる
head(score)## [1] 6 2 6 5 6 4次に、10段階スコアにした官能試験結果のヒストグラムを描きます。
hist(score)
分布がうまく可視化されたように見えるが、「落とし穴」があります。
変数scoreは、10段階（0〜9）のスコアで表現されています。そこで、9段階の各レベルにどれだけのアクセッションが含まれているかをまとめてみます。
t <- table(score)
t## score
## 0 1 2 3 4 5 6 7 8 9
## 2 7 17 13 17 22 18 22 9 6描いたヒストグラムは、クラス全体をうまく表現できていないことがわかります。
上記のように「表」機能を使ってまとめたデータから、棒グラフを描くことができます。
plot(t, xlab = "Liking score", ylab = "Frequency")
棒グラフは「barplot」関数で描くこともできます。
barplot(t, xlab = "Liking score", ylab = "Frequency")
円グラフを描いてみます。
pie(t)
ここからは、変数間の関係を見ていきましょう。
plot(SSC, Overall.Liking)
回帰分析により、データに直線を当てはめる。
plot(SSC, Overall.Liking)
abline(lm(Overall.Liking ~ SSC), col = "green")
このブルーベリーのデータは、5つの品種が3年にわたり、3つの場所で3回ずつ収穫されたデータです。
Genotypeという変数は、品種の違いを表します。品種間で糖度と官能検査の結果の間にどのような関係があるのか、可視化して確認してみましょう。
plot(SSC, Overall.Liking, col = Genotype)
legend("bottomright", levels(Genotype), col = 1:nlevels(Genotype), pch = 1)
ほとんどの品種で、糖度があがると官能試験における評価が高くなることがわかりますが、例外もありそうです（Emerald）。
次に、品種の違いによる糖度や官能試験の評価の違いを箱ひげ図で示してみましょう。
boxplot(SSC ~ Genotype, border = 1:nlevels(Genotype))
boxplot(Overall.Liking ~ Genotype, border = 1:nlevels(Genotype))
SSCが低くても、官能試験の評価が高い品種（Meadowlark）もあり、両者の関係はそれほど単純ではないことが示唆されます。
次に、官能試験の評価と、糖度、酸度のすべての組み合わせについて散布図に描いてみましょう。
x <- data.frame(Overall.Liking, SSC, TA)
pairs(x, col = Genotype)
グラフをファイルに出力する
論文やプレゼンテーションでグラフを使用するために、グラフをPDFファイルに出力できると便利です。ここでは、その方法について簡単に説明します。
先ほど描いた図をfig00.pdfというファイルに出力してみます。
pdf(here("out", "fig00.pdf"))
pairs(x, col = Genotype)
dev.off()## quartz_off_screen
## 2上記のコマンドを実行すると、fig00.pdfというファイルが”out”と名付けられたフォルダ（ディレクトリ）に出力されます。
pdf関数は出力されるグラフの大きさを指定することができます。例えば、pairwise plotを行う変数の数を増やしたいときに、有効です。
pdf(here("out", "fig00_large.pdf"), width = 15, height = 15)
xx <- data.frame(Overall.Liking, Sweetness, Sourness, SSC, TA)
pairs(xx, col = Genotype)
dev.off()## quartz_off_screen
## 2この結果をみると、官能評価における総合評価Overall.Likingと甘さに対する評価Sweetnessに強い相関関係があることが分かります。一方、SweetnessとSSCの相関関係は、Overall.LikingとSweetnessの関係ほどは強くなく、パネリストが感じている甘さは、糖度だけでは十分に説明できていない可能性が示唆されます。
また、複数の図を同じpdfファイルに繰り返し出力すると、複数のページを持つpdfファイルに保存されます。似たような図や関連した図を出力する場合は、1つのpdfファイルにまとめると便利です。
pdf(here("out", "fig00_multi.pdf"))
x <- data.frame(Overall.Liking, SSC, TA)
pairs(x, col = Genotype)
plot(SSC, Overall.Liking, col = Genotype)
abline(lm(Overall.Liking ~ SSC), col = "green")
plot(TA, Overall.Liking, col = Genotype)
abline(lm(Overall.Liking ~ TA), col = "green")
dev.off()## quartz_off_screen
## 2インタラクティブな図を描く
plotlyパッケージを使うと、インタラクティブな図を描くことができます。ここでは、先に描いた図をplotlyの関数plot_lyを用いて描いてみましょう。
まずは、糖度について、ヒストグラムを描きます。
plot_ly(x = SSC, type = "histogram") # plot_ly関数でプロット横向きのヒストグラムも簡単に描くことができます。
plot_ly(y = SSC, type = "histogram")箱ひげ図を描きます。
plot_ly(x = "Soluble Sugar Content", y = SSC, type = "box")品種ごとに箱ひげ図を描きます。
plot_ly(y = SSC, color = Genotype, type = "box")10段階にした官能試験結果のスコアについて、棒グラフを描きます。棒グラフを描く前に、先程と同様に、table関数で頻度を集計しておきます。
t <- table(score)
plot_ly(x = names(t), y = t, type = "bar")円グラフを描きます。
plot_ly(labels = names(t), values = t, type = "pie")次に、2変量（官能試験の評価と糖度）の散布図を描きます。
plot_ly(x = Overall.Liking, y = SSC,
type = "scatter", mode = "markers")品種で色を付けするとともに、ポップアップで示される文字列（Hoverと呼ばれる）に試験地名（Location）を加える。これにより、試験地がどのような影響を与えていそうか確認できる。
plot_ly(x = SSC, y = Overall.Liking,
color = Genotype, text = Location,
type = "scatter", mode = "markers")以下のように、data.frameを使って描くこともできます。
df <- data.frame(Genotype, Location,
Overall.Liking, SSC, TA)
plot_ly(data = df, x = ~TA, y = ~Overall.Liking,
color = ~Genotype, text = ~Location,
type = "scatter", mode = "markers")最後に、3変量間の関係を確認する。
plot_ly(data = df, x = ~SSC, y = ~TA, z = ~Overall.Liking,
color = ~Genotype, text = ~Location,
type = "scatter3d", mode = "markers")確認クイズ
甘さに関連する官能試験評価Sweetnessと、糖度SSC、および、糖成分含量Fructose, Glocose, Sucroseと総含量Total.Sugar間の関係をPDFにpairwise plotとして出力し、これらの変数間の関係を確認してください。
# 自分で書いてみましょう。