# FLEAT VII ワークショップ「外国語教育研究におけるRを用いた統計処理入門」Rコード集

#################### はじめに ####################

# このファイルは、本ワークショップで扱うすべての統計処理を実行するコードをまとめたものです。
# 参考のため載せたコードなど、本ワークショップで扱わないものも含まれますのでご注意ください。
# コードの入力に時間がかかる場合は、このファイルの該当する範囲を選択し、
# コンソールにコピー＆貼り付けしてEnterを押すか、
# 以下のショートカットキーを押してください。コードが実行されます：
# （Rの場合）F5
# （RStudioの場合）Ctrl + Enter



#################### Rの基本 ####################

# 四則演算
3+5
10-3
2*3
100/20

# hakoという名の変数に 5 という数字を代入する
hako <- 5

# hakoという名の変数の中身を確認
hako

# hakoという変数に1〜5を代入
hako <- c(1, 2, 3, 4, 5)
hako

# （参考）代入3パターン
hako <- c(1, 2, 3, 4, 5) 
c(1, 2, 3, 4, 5) -> hako
hako = c(1, 2, 3, 4, 5)



#################### データの読み込み ####################

# 作業ディレクトリの確認
getwd()

# 作業ディレクトリ内の"test.csv"ファイルを読み込み（1行目を見出し行として）
read.csv("test.csv", header= T)	#ファイル名は拡張子（.csv）まで含める

# 読み込んだデータを、datという変数に代入する
dat <- read.csv("test.csv", header = T)

# 変数の中身の確認
dat # すべて表示される
head(dat) # 先頭行から5行目までが表示される

# データフレームから要素を取り出す
# 変数datから事前テストのデータを取り出す場合
dat$pre



#################### パッケージのインストール・読み込み ####################


# パッケージのインストール
install.packages("psych")

# （参考）複数のパッケージを1つの命令でインストールする場合
install.packages(c("psych", "beeswarm"))

# パッケージの読み込み
library(psych)	#パッケージの読み込み



#################### 記述統計量の算出 ####################

# 記述統計量の算出（dat[,1]を引数とする場合）
length(dat[,1])	# データの個数を求める
mean(dat[,1])	# 平均値を求める
median(dat[,1])	# 中央値を求める
max(dat[,1])	# 最大値を求める
min(dat[,1])	# 最小値を求める
var(dat[,1])	# 不偏分散を求める
sd(dat[,1])	# 不偏分散に基づく標準偏差を求める


# summary関数で平均値・第1四分位・中央値・第3四分位・最小値を算出
summary(dat)


# e1071パッケージのskewness関数・kurtosis関数を用いて、歪度・尖度を算出
library(e1071)
skewness(dat[,1])
kurtosis(dat[,1])
# psychパッケージのskew関数・kurt関数を用いて、歪度・尖度を算出
library(psych)
skew(dat[,1])
kurtosi(dat[,1])

# psychパッケージのdescribe関数を用いて記述統計を一通り算出
library(psych)
describe(dat)


# （参考）度数分布表の作成
table(dat[,1]) #身長の度数分布表
table(dat[,2]) #体重の度数分布表


# hist関数でヒストグラム（度数分布図）を描画
hist(dat[,1]) #事前テストのヒストグラム
hist(dat[,2]) #事後テストのヒストグラム


# boxplot関数で箱ひげ図を作成
boxplot(dat)

# 蜂群図（beeswarm plot）を作成
install.packages("beeswarm")
library(beeswarm)
beeswarm(dat)

# 蜂群図を箱ひげ図に重ね描き
boxplot(dat) # 既に箱ひげ図が作成されている状態に
beeswarm(dat, add = T)  # addは重ね描きをするオプション



#################### t検定 ####################

t.test(dat$pre, dat$post, paired = T)
# 1つ目、2つ目の引数にはデータを入れる
# pairedオプションは、対応ありの場合T、対応なしの場合Fを選ぶ

# 効果量（Cohen's d）を求める
install.packages("effsize")
library(effsize)
cohen.d(dat$pre, dat$post)



#################### 一元配置分散分析 ####################

#ANOVA君の読み込み
source("anovakun_483.txt")
source("anovakun_483.txt", encoding = 'CP932')  # Mac/Linuxの方はこちら

# データ読み込み
oneway <- read.csv("oneway.csv", header = T)

# 変数の中身の確認
head(oneway)


anovakun(oneway, "sA", 3, holm = T, eta = T)
# 1つ目の引数はデータ、2つ目の引数は要因計画、3つ目の引数は要因の水準数
# データの並べ方や要因計画の表記、その他オプションについては、以下の井関先生のページを参考にする
# http://riseki.php.xdomain.jp/index.php?ANOVA%E5%90%9B%2FANOVA%E5%90%9B%E3%81%AE%E4%BD%BF%E3%81%84%E6%96%B9
# holmオプションをTにすると、多重比較の方法がHolm法に
# etaオプションをTにすると、効果量（イータ二乗）を算出



#################### 二元配置分散分析 ####################

# データ読み込み
twoway <- read.csv("twoway.csv", header = T)

# 変数の中身の確認
head(twoway)

anovakun(twoway, "AsB", 2, 2, holm = T, peta = T)
# 1つ目の引数はデータ、2つ目の引数は要因計画、3, 4つ目の引数は要因の水準数
# データの並べ方や要因計画の表記、その他オプションについては、以下の井関先生のページを参考にする
# http://riseki.php.xdomain.jp/index.php?ANOVA%E5%90%9B%2FANOVA%E5%90%9B%E3%81%AE%E4%BD%BF%E3%81%84%E6%96%B9
# holmオプションをTにすると、多重比較の方法がHolm法に
# petaオプションをTにすると、効果量（偏イータ二乗）を算出


# （参考）ロング型のデータで二元配置分散分析
# データ読み込み
twoway2 <- read.csv("twoway2.csv", header = T)

# 変数の中身の確認
head(twoway2)

anovakun(twoway2, "AsB", holm = T, peta = T, long = T)
# 1つ目の引数はデータ、2つ目の引数は要因計画
# 水準数は省略可
# データの並べ方や要因計画の表記、その他オプションについては、以下の井関先生のページを参考にする
# http://riseki.php.xdomain.jp/index.php?ANOVA%E5%90%9B%2FANOVA%E5%90%9B%E3%81%AE%E4%BD%BF%E3%81%84%E6%96%B9
# holmオプションをTにすると、多重比較の方法がHolm法に
# petaオプションをTにすると、効果量（偏イータ二乗）を算出



#################### 相関分析 ####################

# データ読み込み
soukan <- read.csv("soukan.csv", header=T)	#データの読み込み

# 変数の中身の確認
head(soukan)

# 散布図
plot(soukan)

# 相関係数（行列）の算出
cor(soukan)

# 相関係数・無相関検定・相関係数の信頼区間を算出する
cor.test(soukan$read, soukan$vocab, method = "pearson")
# methodが"spearman", "kendall"の場合は、信頼区間は算出されない

# （参考）スピアマンの順位相関係数の場合は、以下のようにmethod = "pearson"と、rank関数を組み合わせると、信頼区間が算出される(表示はピアソンのまま)
cor.test(rank(soukan$read), rank(soukan$vocab), method = "pearson")



#################### 回帰分析 ####################

# データ読み込み
kaiki <- read.csv("kaiki.csv", header = T)

# 変数の中身の確認
head(kaiki)

# 相関行列の確認
cor(kaiki)

# 単回帰分析
res1 <- lm(kaiki$toeic ~ kaiki$read)
summary(res1)

# 重回帰分析
res2 <- lm(kaiki$toeic ~ kaiki$read + kaiki$listen + kaiki$vocab)
summary(res2)

# 多重共線性の指標 VIFを算出（2未満だと問題ないといわれる） 
library(car)
vif(res2)