生物統計学2019 第三回 授業におけるＲコード
Gakushi Ishimura
October 29, 2019
生物統計学入門 2019 石村 学志
もし、Ｒをインストールしていない場合は https://ja.wikipedia.org/wiki/R言語 を読んで自分のパソコンにインストールしてください。 統計科学研究所のサイト http://www.statistics.co.jp/reference/software_R/free_software-R.htm に使い方と簡単な使用方法が掲載されています。また、Rstudioもインストールすることを推奨します。https://qiita.com/hujuu/items/ddd66ae8e6f3f989f2c0 にRの解説とともに、Rstudioのインストールについて書いてあります。
分析に必要なライブラリーのインストール
R は関数とデータを機能別に分類してパッケージ（library）として容易されており、 必要に応じてダウンロードなどで読み込む必要がある。
library(dplyr) # data scienceでの新しい概念である"tidy data"のためのパッケージ
library(DT)フィールド用データ
使用するサンプルデータサイトは’R for Fisheries Analyses’です。 http://derekogle.com/fishR/ http://derekogle.com/fishR/data/byTopic
‘R’をインストールしたあとで ’FSdata’‘FAS’’dplyr’をインストールしてください パッケージのインストール方法は http://kawamurakaeru.lolipop.jp/scw/blog/2016/02/05/rstudio でパッケージインストールを参照してください。
FSAとFSAdataは水産資源分析のパッケージとデータです。 dplyrは、整然データ(tidy data)という新しいデータサイエンス分野における取り扱いやすいデータの形を つくるためのパッケージです。 詳しくは下記参照 http://id.fnshr.info/2017/01/09/tidy-data-intro/
library(dplyr) # data scienceでの新しい概念である"tidy data"のためのパッケージ
library(FSA) # 水産資源分析パッケージ
library(FSAdata) # FSAのためのデータフィールド系データ：チリのChinook Salmon（キングサーモン）データ ”Lengths and weights for Chinook Salmon from three locations in Argentina” https://www.rdocumentation.org/packages/FSA/versions/0.8.16/topics/ChinookArg
実験系データ：モルモットの2種類のビタミン投与による歯の成長データ “The response is the length of odontoblasts (cells responsible for tooth growth) in 60 guinea pigs. Each animal received one of three dose levels of vitamin C (0.5, 1, and 2 mg/day) by one of two delivery methods, orange juice or ascorbic acid (a form of vitamin C and coded as VC)” https://www.rdocumentation.org/packages/datasets/versions/3.5.1/topics/ToothGrowth
# どのようなデータか説明をみる
help.search("ChinookArg",package=c("FSAdata","FSA"))
help.search("ToothGrowth",package=c("datasets"))
# データを加工するためにWorking_Dataという変数にデータを格納する。
Working_Data=ChinookArg
Working_Data_2=ToothGrowth
## 変数名の表示
names(Working_Data)## [1] "tl" "w" "loc"names(Working_Data_2)## [1] "len" "supp" "dose"## データの大きさ
dim(Working_Data)## [1] 112 3dim(Working_Data_2)## [1] 60 3## データ全体の概観をみる
glimpse(Working_Data)## Observations: 112
## Variables: 3
## $ tl <dbl> 120.1, 115.0, 111.2, 110.2, 110.0, 109.7, 105.0, 100.1, 98.…
## $ w <dbl> 17.9, 17.2, 16.8, 15.8, 14.3, 13.8, 12.8, 11.7, 12.8, 14.8,…
## $ loc <fct> Argentina, Argentina, Argentina, Argentina, Argentina, Arge…glimpse(Working_Data_2)## Observations: 60
## Variables: 3
## $ len <dbl> 4.2, 11.5, 7.3, 5.8, 6.4, 10.0, 11.2, 11.2, 5.2, 7.0, 16.5…
## $ supp <fct> VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC…
## $ dose <dbl> 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 1.0…### または
str(Working_Data)## 'data.frame': 112 obs. of 3 variables:
## $ tl : num 120 115 111 110 110 ...
## $ w : num 17.9 17.2 16.8 15.8 14.3 13.8 12.8 11.7 12.8 14.8 ...
## $ loc: Factor w/ 3 levels "Argentina","Petrohue",..: 1 1 1 1 1 1 1 1 1 1 ...str(Working_Data_2)## 'data.frame': 60 obs. of 3 variables:
## $ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
## $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
## $ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...datatable(ChinookArg)glimpse(ChinookArg)## Observations: 112
## Variables: 3
## $ tl <dbl> 120.1, 115.0, 111.2, 110.2, 110.0, 109.7, 105.0, 100.1, 98.…
## $ w <dbl> 17.9, 17.2, 16.8, 15.8, 14.3, 13.8, 12.8, 11.7, 12.8, 14.8,…
## $ loc <fct> Argentina, Argentina, Argentina, Argentina, Argentina, Arge…次に使いやすいように変数名を書き換える。
names(Working_Data)=c("Length","Weight","Location")
## 変数名の表示で名前が変わっているかチェックする。
names(Working_Data)## [1] "Length" "Weight" "Location"names(Working_Data_2)=c("ToothLength","Supplement","Dose")
## 変数名の表示で名前が変わっているかチェックする。
names(Working_Data_2)## [1] "ToothLength" "Supplement" "Dose"３つの場所の名前を確認
unique(Working_Data$Location)## [1] Argentina Petrohue Puyehue
## Levels: Argentina Petrohue Puyehueとりあえず最初の2つと最後の5つのデータをみてみる。
head(Working_Data,n=2)## Length Weight Location
## 1 120.1 17.9 Argentina
## 2 115.0 17.2 Argentinatail(Working_Data,n=5)## Length Weight Location
## 108 32.1 2.8 Puyehue
## 109 31.9 0.3 Puyehue
## 110 29.2 0.3 Puyehue
## 111 25.2 0.3 Puyehue
## 112 18.0 0.1 Puyehueデータの叙述統計(discriptive statistics)をみる。
summary(Working_Data)## Length Weight Location
## Min. : 18.00 Min. : 0.100 Argentina:34
## 1st Qu.: 65.60 1st Qu.: 2.975 Petrohue :30
## Median : 81.60 Median : 6.400 Puyehue :48
## Mean : 78.86 Mean : 6.678
## 3rd Qu.: 94.67 3rd Qu.: 9.475
## Max. :120.10 Max. :17.900体長のデータをヒストグラムでみてみると。
hist(Working_Data$Length)
体重のデータをヒストグラムでみてみると。
hist(Working_Data$Weight)
簡単にするためにデータのなかかから場所・アルゼンチンだけでのデータをとりだす。
Working_Data2=filter(Working_Data,Working_Data$Location=='Argentina')データの叙述統計(discriptive statistics)をみる。
summary(Working_Data2)## Length Weight Location
## Min. : 59.90 Min. : 3.900 Argentina:34
## 1st Qu.: 82.90 1st Qu.: 7.425 Petrohue : 0
## Median : 92.05 Median : 9.400 Puyehue : 0
## Mean : 91.31 Mean :10.200
## 3rd Qu.: 98.75 3rd Qu.:12.575
## Max. :120.10 Max. :17.900体長のデータをヒストグラムでみてみると。
hist(Working_Data2$Length)
体重のデータをヒストグラムでみてみると。
hist(Working_Data2$Weight)
体長と伸張のデータの関係をプロットしてみると。
plot(Working_Data2$Length,Working_Data2$Weight)
体長と体重のデータの関係をもうちょっと綺麗にプロットしてみると。
par(family = "HiraKakuProN-W3") # Rの出力に日本語表示するための環境設定
plot(Working_Data2$Length,Working_Data2$Weight,
xlim=c(0,max(Working_Data$Length)*1.2),
ylim=c(0,max(Working_Data$Weight)*1.2),
xlab="体長",ylab="体重")
体長と体重の関係をみる
output.lm <- lm(Working_Data2$Length~Working_Data2$Weight)
summary(output.lm)##
## Call:
## lm(formula = Working_Data2$Length ~ Working_Data2$Weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.942 -4.835 -1.391 5.126 15.880
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58.6394 3.8580 15.199 3.43e-16 ***
## Working_Data2$Weight 3.2029 0.3562 8.992 2.85e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.571 on 32 degrees of freedom
## Multiple R-squared: 0.7165, Adjusted R-squared: 0.7076
## F-statistic: 80.86 on 1 and 32 DF, p-value: 2.851e-10