生物統計学入門　2020 石村　学志

# libraryの読み込み
library(tidyverse)
## ─ Attaching packages ──────────────────── tidyverse 1.3.0 ─
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.0
## ─ Conflicts ───────────────────── tidyverse_conflicts() ─
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(readr)　　#　csv file の読み込み
library(readxl)　#  xls file　の読み込み

# "R"が今、参照している場所(ディレクトリーもしくはフォルダー)の確認
getwd()
## [1] "/Users/gakugaku/Desktop/BioStat2020"
# "R"が参照してい場所（ディレクトリー）を変更
setwd("/Users/gakugaku/Desktop/BioStat2020/")


#　参照しているディレクトリーにあるファイルをみる
dir()
##  [1] "BioStat2020"                 "BioStat2020_1.html"         
##  [3] "BioStat2020_1.Rmd"           "BioStat2020_Lec4_cache"     
##  [5] "BioStat2020_Lec4.html"       "BioStat2020_Lec4.log"       
##  [7] "BioStat2020_Lec4.Rmd"        "BioStat2020_Lec4.tex"       
##  [9] "ChinookDataBioStat2020.csv"  "ChinookDataBioStat2020.xlsx"
## [11] "Lec4.html"                   "Lec4.Rmd"                   
## [13] "Untitled_files"              "Untitled.html"              
## [15] "Untitled.Rmd"
#   csvファイルの読み込み
ChinookData_CSV <- read_csv("ChinookDataBioStat2020.csv")
## 
## ─ Column specification ────────────────────────────
## cols(
##   Length = col_double(),
##   Weight = col_double(),
##   Location = col_character()
## )
#   エクセルファイルの読み込み
ChinookData_XLS <- read_excel("ChinookDataBioStat2020.xlsx")


#   Working_Data という違うデータセットを作る
Working_Data <- ChinookData_CSV 

#  変数名の表示
names(Working_Data)
## [1] "Length"   "Weight"   "Location"
## データの大きさ
dim(Working_Data)
## [1] 112   3
##  データ全体の概観をみる
glimpse(Working_Data)
## Rows: 112
## Columns: 3
## $ Length   <dbl> 120.1, 115.0, 113.8, 112.9, 111.2, 110.2, 110.0, 109.7, 109…
## $ Weight   <dbl> 17.9, 17.2, 15.0, 16.0, 16.8, 15.8, 14.3, 13.8, 11.3, 13.3,…
## $ Location <chr> "Kamaishi", "Kamaishi", "Miyako", "Miyako", "Kamaishi", "Ka…
###  または
str(Working_Data)
## tibble [112 × 3] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ Length  : num [1:112] 120 115 114 113 111 ...
##  $ Weight  : num [1:112] 17.9 17.2 15 16 16.8 15.8 14.3 13.8 11.3 13.3 ...
##  $ Location: chr [1:112] "Kamaishi" "Kamaishi" "Miyako" "Miyako" ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   Length = col_double(),
##   ..   Weight = col_double(),
##   ..   Location = col_character()
##   .. )

使用するサンプルデータサイトは’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

#  どのようなデータか説明をみる
help.search("ChinookArg",package=c("FSAdata","FSA"))

#  データを加工するためにWorking_Dataという変数にデータを格納する。
Working_Data <- ChinookArg  

##  変数名の表示
names(Working_Data)
## [1] "tl"  "w"   "loc"
##  データの大きさ
dim(Working_Data)
## [1] 112   3
##  データ全体の概観をみる
glimpse(Working_Data)
## Rows: 112
## Columns: 3
## $ tl  <dbl> 120.1, 115.0, 111.2, 110.2, 110.0, 109.7, 105.0, 100.1, 98.0, 92…
## $ w   <dbl> 17.9, 17.2, 16.8, 15.8, 14.3, 13.8, 12.8, 11.7, 12.8, 14.8, 9.7,…
## $ loc <fct> Argentina, Argentina, Argentina, Argentina, Argentina, Argentina…
###  または
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 ...

次に使いやすいように変数名を書き換える。

Working_Data <- ChinookData_CSV 
names(Working_Data) <- c("Length","Weight","Location")
##  変数名の表示で名前が変わっているかチェックする。
names(Working_Data)
## [1] "Length"   "Weight"   "Location"
# データの型を個別に確認
class(Working_Data$Length)
## [1] "numeric"
class(Working_Data$Weight)
## [1] "numeric"
class(Working_Data$Location)
## [1] "character"

３つの場所の名前を確認

unique(Working_Data$Location)
## [1] "Kamaishi" "Miyako"   "Rikuzen"

とりあえず最初の2つと最後の5つのデータをみてみる。

head(Working_Data,n=2)
tail(Working_Data,n=5)

データの叙述統計(discriptive statistics)をみる。

summary(Working_Data)
##      Length           Weight         Location        
##  Min.   : 18.00   Min.   : 0.100   Length:112        
##  1st Qu.: 65.60   1st Qu.: 2.975   Class :character  
##  Median : 81.60   Median : 6.400   Mode  :character  
##  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　==　'Kamaishi')

データの叙述統計(discriptive statistics)をみる。

summary(Working_Data2)
##      Length           Weight         Location        
##  Min.   : 59.90   Min.   : 3.900   Length:34         
##  1st Qu.: 82.90   1st Qu.: 7.425   Class :character  
##  Median : 92.05   Median : 9.400   Mode  :character  
##  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)

体長と体重のデータの関係をもうちょっと綺麗にプロットしてみると。

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="体重")

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

より美しい絵を描くためにggplotを使う

###  ヒストグラムを描く (数値データ(numeric) のみ)
g1　<- ggplot(Working_Data, aes(x = Length,fill=Location)) 
g2  <- g1+geom_histogram(alpha = 0.5,position="stack")
plot(g2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

g3　<- ggplot(Working_Data, aes(x = Length,fill=Location)) 
g4  <- g3+geom_histogram(alpha = 0.5,position="identity")
plot(g4)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

###  体長と体重の関係を描く
g5 <- ggplot(Working_Data, aes(x = Length,y=Weight, col = Location)) 
g6 <- g5 +geom_point()  
plot(g6)