目次1

目次2

目次3

このまま入力すれば本文となる。
改行する時は文末に半角ダブルスペースを加えるか、一行以上あける。

イタリック
太字　　

リストは
このように
入力する

アルファベットの
リスト

数字の
リスト

LaTeXでの数式入力は$Y= \beta_0+\beta_1x_1+\beta_2x_2$と入力する。

コードを本文中に記載する場合はcodeとする。

Rmdファイルを見やすいフォーマットのファイルに変換することをKnit（ニット）と言う。基本的にknitしたrmdファイルは自動で上書き保存される。rmdファイル内にエラーを生み出すようなコードがある場合、knitはうまくいかない。該当コードをコメントアウトするか、後に述べるerror statementを使えば回避できる。

Rmdファイルのメリットは、コードだけでなくコードを走らせた後の結果や図も表示されることである。チーム間で解析をシェアする時に何をしているか分かりやすい、再現しやすいという強みがある。

# コードはこのようなチャンクの中で書く
# でコメントアウトできる

{r}の中でさまざまな設定ができる。
(例)
{r チャンク名} チャンクに名前をつけることでデバッグがしやすくなる。
echo=T でチャンク内のコードを表示、echo=Fで非表示にできる。　
result=F で結果を非表示にできる。　　
error=T でエラーが起きたときもknitできるようにする。
message=F でメッセージを非表示にできる。
warning=F でwarningを非表示にできる。
eval=F はコードの表示はするが、コードを走らせることはしない。
一括で設定を変えたい（グローバル設定をする）ときはknitr::opts_chunk$set(echo = TRUE)の中に上記の条件を入力する。

Rコードの基礎

解析で使う頻度の高いコードについて簡単に紹介する。

Package/library

Rに元々実装されている機能のことをbase Rと呼ぶ。データの取り込みや解析を行う際にはbaseでは足りないことが多く、パッケージというものを呼び起こす必要がある。調べたコードを実行する時、しばしば関数が見つからないというエラーが出る。その際は”使いたい関数 + r + package”で調べるとインストールすべきパッケージがわかる。

# パッケージがインストールされていない時
# install.packages("tidyverse")
# install.packages("ggplot2")

# パッケージの呼び出し
library(tidyverse) # data manipulation 
library(ggplot2) # data visualization
library(readr) # import csv and tsv
library(readxl) # import excel

ちなみにrmdに画像を貼ることもできる。

基本演算

変数名を決定する際のルール、慣習：

数字からはじめない
スペースを設けない (blood pressureではなくblood_pressureとする)
_, -, .以外の記号(@, /, *, &, $ など)は避ける
大文字、小文字を統一する (小文字とすることが多い)
できるだけ短く、informativeなものを

# 変数の定義
weight <- 50
height <- 1.67

# 演算
bmi <- round(weight / height^2, 2) # 小数第2位まで表示
print(bmi)
## [1] 17.93

# 電卓のように数字をそのまま打ち込んでもいい
1+2
## [1] 3
100-7
## [1] 93
2*3
## [1] 6
10/5
## [1] 2
10^2 # 指数の計算
## [1] 100
exp(0) # ネイビア数の計算
## [1] 1
log(exp(1)) # 対数の計算 (指定しなければbaseはネイピア数)
## [1] 1

リストとデータフレーム

# リスト
list1 <- c("a", "b", "c")
list2 <- c(1, 2, 3)
list3 <- c(2, 2, 2)

# リストの計算
list2 + list3
## [1] 3 4 5
list2 * list3
## [1] 2 4 6

# データフレーム
df <- data.frame(
  "a" = list2, # 左辺が列名で右辺が対応する列の値
  "b" = list3
)
print(df)
##   a b
## 1 1 2
## 2 2 2
## 3 3 2

# 変数タイプの確認
class(list1)
## [1] "character"
class(list2)
## [1] "numeric"
class(df)
## [1] "data.frame"
class(TRUE)
## [1] "logical"

# 変数タイプを無理矢理指定することもできる
as.character(list2)
## [1] "1" "2" "3"
class(as.character(list2))
## [1] "character"

データの取り込み

xlsx: df <- read_xlsx(file = "PATH")
csv: df <- read_csv(file = "PATH", col_names = TRUE)

データハンドリング (`dplyr`)

RやPythonなどの言語の最大の強みは、元データを保ちながらデータを編集、解析し、その記録がすべて残ることである（reproducibility)。

df <- iris

# データサマリー
head(iris) # 最初の数行を確認
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
summary(iris) # meanやmedianなどをざっくり確認
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
## 

# データの呼び起こし
colnames(df) # 列名のリストを表示
## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width"  "Species"
df$Sepal.Length |>  # 特定の列を指定
  head()
## [1] 5.1 4.9 4.7 4.6 5.0 5.4
length(df$Sepal.Length) # 行の総数
## [1] 150
length(df) # 列の総数
## [1] 5


# 変数名を整える便利な機能
df |> 
  janitor::clean_names() |> 
  head()
##   sepal_length sepal_width petal_length petal_width species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

# その他のデータハンドリング
df |> 
  select(c(Sepal.Length, Sepal.Width)) |> # 特定の列をデータフレームから抽出する 
  filter(Sepal.Length <= 5) |> # Sepal.Lengthが5以下となる行のみ抽出する
  mutate(
    Sepal.Length = c(1:length(Sepal.Length)), # 特定の行の値を書き換える
    additional = rnorm(n = length(Sepal.Length), mean = 0, sd = 5) # 新しい変数を追加する（rnormは正規分布に従ったランダムな値を生成する）
  ) |> 
  rename(
    "width" = Sepal.Width # 変数名を変更する
  ) |> 
  head()
##   Sepal.Length width additional
## 1            1   3.0  4.3749603
## 2            2   3.2 -5.2128220
## 3            3   3.1 -1.3092337
## 4            4   3.6 -0.3410321
## 5            5   3.4 -4.4258071
## 6            6   3.4 -3.4825358

Visualization (`ggplot2`)

df <- iris

# ヒストグラム（連続変数）
ggplot(df, aes(x = Petal.Length)) +
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


# 棒グラフ（カテゴリカル変数）
ggplot(df, aes(x = Species)) +
  geom_bar()


# 箱ひげ図
ggplot(df, aes(x = Species, y = Petal.Length)) + # x軸はカテゴリカル変数
  geom_boxplot()


# 分散図
ggplot(df, aes(x = Petal.Width, y = Petal.Length)) + 
  geom_point(aes(color = Species), alpha = .5) + # colorによって第三のカテゴリカル変数を図示できる
  labs( # 図中のラベル表記を指定する
    title = "Scatter Plot",
    x = "Petal Width",
    y = "Petal Length",
    color = "Species",
    caption = "Data from iris"
  )

基本的な解析

カイ2乗検定

Sepal.Lengthをbinary（High/Low）に分類し、Speciesの違いと関連があるかを調べる。

# データ準備
# categorize Sepal.Length into "High" and "Low" based on the median
iris$Sepal.Length.category <- ifelse(iris$Sepal.Length > median(iris$Sepal.Length), "High", "Low") # ifelse(A > B, C, D): もしA>Bが真であればCとし、偽であればDとする
table(iris$Sepal.Length.category, iris$Species) # 各SpeciesにHighとLowがどれだけあるか集計して表にする
##       
##        setosa versicolor virginica
##   High      0         26        44
##   Low      50         24         6

# カイ2乗検定でcategoryとSpeciesの関係を調べる
chisq.test(table(iris$Sepal.Length.category, iris$Species)) |> 
  broom::tidy() |> # 結果が表としてリターンされる (必須ではない)
  knitr::kable(digits = 2) # knitするときに結果がきちんと表として表示される。digitsで小数何桁まで表示するか選べる

statistic	p.value	parameter	method
78.64	0	2	Pearson’s Chi-squared test

T検定

2種類のSpecies間でSepal.Lengthの平均の差が統計的に優位か調べてみる。

# データ準備
df <- iris
# subset data for two species
subset_iris <- subset(iris, Species %in% c("setosa", "versicolor")) # subsetはselectと同義。 %in%は含まれているかどうか調べる時に使う

# perform t-test on Sepal.Length between setosa and versicolor
t.test(Sepal.Length ~ Species, data = subset_iris) |> 
  broom::tidy() |>
  knitr::kable(digits = 2)

estimate	estimate1	estimate2	statistic	p.value	parameter	conf.low	conf.high	method	alternative
-0.93	5.01	5.94	-10.52	0	86.54	-1.11	-0.75	Welch Two Sample t-test	two.sided

ANOVA

3種類のSpecies間でSepal.Lengthの平均の差が統計的に優位か調べてみる。

# データ準備
df <- iris

# perform ANOVA to test if Sepal.Length differs across Species
aov(Sepal.Length ~ Species, data = iris) |> 
  broom::tidy() |>
  knitr::kable(digits = 2)

term	df	sumsq	meansq	statistic	p.value
Species	2	63.21	31.61	119.26	0
Residuals	147	38.96	0.27	NA	NA

Linear Regression

Sepal.Lengthを結果、Sepal.Widthを予測変数としてモデルを組んでみる。

# linear regression predicting Sepal.Length using Sepal.Width
linear_model <- lm(Sepal.Length ~ Sepal.Width, data = iris)
summary(linear_model)
## 
## Call:
## lm(formula = Sepal.Length ~ Sepal.Width, data = iris)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5561 -0.6333 -0.1120  0.5579  2.2226 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.5262     0.4789   13.63   <2e-16 ***
## Sepal.Width  -0.2234     0.1551   -1.44    0.152    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8251 on 148 degrees of freedom
## Multiple R-squared:  0.01382,    Adjusted R-squared:  0.007159 
## F-statistic: 2.074 on 1 and 148 DF,  p-value: 0.1519

# またはこれまでと同様に下記のように出力しても良い
lm(Sepal.Length ~ Sepal.Width, data = iris) |> 
  broom::tidy() |>
  knitr::kable(digits = 2)

term	estimate	std.error	statistic	p.value
(Intercept)	6.53	0.48	13.63	0.00
Sepal.Width	-0.22	0.16	-1.44	0.15


# 各種結果を呼び出すことも可能
linear_model$coefficients 
## (Intercept) Sepal.Width 
##   6.5262226  -0.2233611

# 正規性を可視化する
qqnorm(linear_model$residuals)
qqline(linear_model$residuals, col = "red", lwd = 2) # reference line for normality

Logistic Regression

setosaとversicolorを結果、Sepal.LengthとSepal.Widthを予測変数としてモデルを組んでみる。

# subset the data to include only two species for binary outcome
iris_binary <- subset(iris, Species %in% c("setosa", "versicolor"))

# convert Species to a binary factor
iris_binary$Species <- factor(iris_binary$Species, levels = c("setosa", "versicolor")) # factorはbinary変数にダミー変数を割り当てる

# git logistic regression model
logit_model <- glm(Species ~ Sepal.Length + Sepal.Width, 
                   data = iris_binary, 
                   family = binomial()) # 二項分布のほかにポワソン分布、正規分布、ガンマ分布などを指定できる

summary(logit_model)
## 
## Call:
## glm(formula = Species ~ Sepal.Length + Sepal.Width, family = binomial(), 
##     data = iris_binary)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -360.6   195972.9  -0.002    0.999
## Sepal.Length    131.8    64577.0   0.002    0.998
## Sepal.Width    -110.1    55361.5  -0.002    0.998
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1.3863e+02  on 99  degrees of freedom
## Residual deviance: 7.1185e-09  on 97  degrees of freedom
## AIC: 6
## 
## Number of Fisher Scoring iterations: 25

# または
glm(Species ~ Sepal.Length + Sepal.Width, 
                   data = iris_binary, 
                   family = binomial()) |> 
  broom::tidy() |>
  knitr::kable(digits = 2)

term	estimate	std.error	p.value
(Intercept)	-360.61	195972.86	1
Sepal.Length	131.79	64576.99	1
Sepal.Width	-110.13	55361.50	1

最後に

Rはフレキシブルな解析ツールであり、日々新しいパッケージが開発されている。最近では機械学習を用いた論文をよく目にするようになったが、それらのほとんどがRやPythonによって解析されている。また、ChatGPTなどとも相性が良く、ラーニングカーブは下がってきている印象。利点が多い一方で、SASやSPSS、STATAと異なり第三者機関の検閲がないため特に新しく開発されたパッケージ(tidymodelsなど)では稀に間違いが指摘されることがあるため、注意されたい。

便利な参考資料

Rによる統計入門
 R for Data Science
R Graphics Cookbook, 2nd edition
Table 1を出力できるgtsummary
チートシート一覧
 RStudioとGitHubとの連携

Rの基礎

Yuki Joyama

2024-10-31

目次1

目次2

目次3

Rコードの基礎

Package/library

基本演算

リストとデータフレーム

データの取り込み

データハンドリング (`dplyr`)

Visualization (`ggplot2`)

基本的な解析

カイ2乗検定

T検定

ANOVA

Linear Regression

Logistic Regression

最後に

便利な参考資料

Rの基礎

Yuki Joyama

2024-10-31

目次1

目次2

目次3

Rコードの基礎

Package/library

基本演算

リストとデータフレーム

データの取り込み

データハンドリング (dplyr)

Visualization (ggplot2)

基本的な解析

カイ2乗検定

T検定

ANOVA

Linear Regression

Logistic Regression

最後に

便利な参考資料

データハンドリング (`dplyr`)

Visualization (`ggplot2`)