大変勉強になりました。

列ごとに、ファイルの値（NAがある）が2SDより高い値をNAにしたいです。

Data01 <-read_excel("raw.xlsx",1)

Data02 <-apply(Data01, 2, sd, na.rm = TRUE)

Data01[Data01 > 2*Data02]<- NA

以上のコードをして最後は列ごとにうまく出来なかったです。

ご提言があれば大変嬉しいです。

Data01 <-read_excel("raw.xlsx",1)

set.seed(2020)
a <- c(rnorm(998,10,1),NA,30)
b <- c(rnorm(998,20,1),NA,40)
c <- c(rnorm(998,30,1),NA,50)
Data01 <- data.frame(a,b,c)
head(Data01)

Data02 <-apply(Data01, 2, sd, na.rm = TRUE)

Data01_mean <-apply(Data01, 2, mean, na.rm = TRUE)
Data01_sd <-apply(Data01, 2, sd, na.rm = TRUE)
Data01_mean
Data01_sd

low <- Data01_mean - 2*Data01_sd
high <- Data01_mean + 2*Data01_sd
low
high

Data01$a[(Data01$a < low["a"]) | (high["a"] < Data01$a)]

Data02 <- Data01
Data02$a[(Data02$a < low["a"]) | (high["a"] < Data02$a)] <- NA
Data02$b[(Data02$b < low["b"]) | (high["b"] < Data02$b)] <- NA
Data02$c[(Data02$c < low["c"]) | (high["c"] < Data02$c)] <- NA

#もし初めての場合はインストール
install.packages("tidyverse")
#インストール済ならlibrary関数で呼び出す
library(tidyverse)

low_high <-
  summarize(Data01,
          a_low = mean(a, na.rm = TRUE)-2*sd(a, na.rm = TRUE),
          a_high = mean(a, na.rm = TRUE)+2*sd(a, na.rm = TRUE),
          b_low = mean(b, na.rm = TRUE)-2*sd(b, na.rm = TRUE),
          b_high = mean(b, na.rm = TRUE)+2*sd(b, na.rm = TRUE),
          c_low = mean(c, na.rm = TRUE)-2*sd(c, na.rm = TRUE),
          c_high = mean(c, na.rm = TRUE)+2*sd(c, na.rm = TRUE)) 
class(low_high) #型を確認するとdata.frameとなっている
low_high <- unlist(low_high)
class(low_high) #型を確認するとnumeric（数値型）となっている
low_high

#%>%使う場合
low_high <-
  summarize(Data01,
          a_low = mean(a, na.rm = TRUE)-2*sd(a, na.rm = TRUE),
          a_high = mean(a, na.rm = TRUE)+2*sd(a, na.rm = TRUE),
          b_low = mean(b, na.rm = TRUE)-2*sd(b, na.rm = TRUE),
          b_high = mean(b, na.rm = TRUE)+2*sd(b, na.rm = TRUE),
          c_low = mean(c, na.rm = TRUE)-2*sd(c, na.rm = TRUE),
          c_high = mean(c, na.rm = TRUE)+2*sd(c, na.rm = TRUE)) %>%
  unlist()
low_high

Data02 <- 
  Data01 %>% 
  mutate(a = ifelse((low_high["a_low"] < a) & (a < low_high["a_high"]), a, NA),
         b = ifelse((low_high["b_low"] < b) & (b < low_high["b_high"]), b, NA),
         c = ifelse((low_high["c_low"] < c) & (c < low_high["c_high"]), c, NA))

drop_2sd <- function(x){
  ifelse((mean(x, na.rm = TRUE) - 2 * sd(x, na.rm = TRUE) < x) & (x < mean(x, na.rm = TRUE) + 2 * sd(x, na.rm = TRUE)), x, NA)
}

Data02 <- apply(Data01, 2, drop_2sd)

第4章は統計を扱います。

今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。

シロート統計学はEZRを使った統計分析をわかりやすく解説されています。

第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか？といった視点で進めていきます。

今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください！

今回はFriedman（フリードマン）検定を紹介します

download.file(url = “ファイルのURL”,

　　　　　　　 destfile = “保存したい場所/ファイル名”)
urlはデモデータで右クリック → リンクのアドレスをコピー
destfileは保存場所と保存のファイル名を指定します。

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-friedman-rank-sum-test.xlsx"
destfile = "data/demo-friedman-rank-sum-test.xlsx"

download.file(url, destfile)

4.データの読み込み

データを読み込みます。
今回は【4-0】第4章を進めていく上での準備で行った方法で進めます。

View Fileでデータを確認します。

library(readxl)
grip_friedman <- read_excel("data/demo-friedman-rank-sum-test.xlsx")
View(grip_friedman)

friedman.test(cbind(grip_friedman$W0,grip_friedman$W1,grip_friedman$W2))

id：30人の内誰のデータかわかるためのid（氏名でも可）
time：W0, W1, W2の3つ
grip：今回測定した握力のデータ

library(tidyverse)

#gather関数の場合
grip_friedman_long <-
  grip_friedman %>% 
  rowid_to_column("id") %>% 
  gather(key = time, value = grip, W0:W2, factor_key = TRUE)

#pivot_longer関数の場合
grip_friedman_long <-
  grip_friedman %>% 
  rowid_to_column("id") %>% 
  pivot_longer(cols = c(W0,W1,W2), names_to = "time", values_to = "grip", names_ptypes = list(time = factor())) 

friedman.test(grip ~ time | id, data = grip_friedman_long)

#Bonferroniの方法
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "bonferroni", paired = TRUE, exact = FALSE)

#Holmの方法
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "holm", paired = TRUE, exact = FALSE)

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-friedman-rank-sum-test.xlsx"
destfile = "data/demo-friedman-rank-sum-test.xlsx"

download.file(url, destfile)

library(readxl)
grip_friedman <- read_excel("data/demo-friedman-rank-sum-test.xlsx")
View(grip_friedman)

friedman.test(cbind(grip_friedman$W0,grip_friedman$W1,grip_friedman$W2))

#pivot_longer関数の場合
grip_friedman_long <-
  grip_friedman %>% 
  rowid_to_column("id") %>% 
  pivot_longer(cols = c(W0,W1,W2), names_to = "time", values_to = "grip", names_ptypes = list(time = factor())) 

#gather関数の場合
grip_friedman_long <-
  grip_friedman %>% 
  rowid_to_column("id") %>% 
  gather(key = time, value = grip, W0:W2, factor_key = TRUE)

friedman.test(grip ~ time | id, data = grip_friedman_long)

pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "bonferroni", paired = TRUE, exact = FALSE)
pairwise.wilcox.test(grip_friedman_long$grip, grip_friedman_long$time, p.adjust.method = "holm", paired = TRUE, exact = FALSE)
 

第4章は統計を扱います。

今回も「シロート統計学」のハルさんとコラボレーションすることとなりました。

シロート統計学はEZRを使った統計分析をわかりやすく解説されています。

第4章はシロート統計学で使われていたEZRをRで行うとどうなるのか？といった視点で進めていきます。

今回使うデータもハルさんのサイトと同じものを使わせでいただく事になりました。それぞれ見比べることで参考にしてみてください！

まずはハルさんのサイトをご覧ください。
ここで行っている解析をRで実施する方法を紹介します。


1.データの準備

今回は前回のデータを使います。記事の1~4までをご準備ください。


もし前回の記事を見るなんて面倒くさい！という方は以下のコードを実行してください。
（このコードはデータの保存場所が今までのの記事と違います）

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-repeated-measures-anova.xlsx"
destfile <- "demo-repeated-measures-anova.xlsx" #今まで第4章の記事を読まれて実践されてる方は destfile <- "data/demo-repeated-measures-anova.xlsx"

download.file(url, destfile)

#データの読み込み
library(readxl)
grip_rep_anova <- read_excel("data/demo-repeated-measures-anova.xlsx", 
                             range = "B2:E32")
View(grip_rep_anova)

前回はある運動プログラムを行い、0週目・1週目・2週目で握力を測定した、ということを想定した仮想データとなっています。今回は握力の経時的検討に加えて男女の群分けも行います。

2.被験者間要因と被験者内要因と分散分析のデザインについて
分散分析を行う上で大切なのがｖについてです。

被験者間要因はいわゆる対応のないデータのことで被験者全体をA法、B法の2つに割り振るといった方法です。なので同じ人がA法とB法を行うことはありません。
被験者内要因はいわゆる対応のあるデータのことで同じ人が繰り返し測定します。

そして被験者間要因は列の左側に、被験者内要因は列の右側に並べると後々分析しやすくなります。
また被験者間要因は縦に、被験者内要因は横につなげるとEZRで分析を行うには都合がいいです。

仮に下のデータがあったとします。

この並びが分散分析を行う上での基本的な形になります。

要因が3つ（sex, treatment, 繰り返し）あるのでABC
被験者間要因と被験者内要因の間にsを付ける

上記の場合はABsCデザインと言います。
また被験者間要因と被験者内要因が両方あるデザインを混合計画とも呼んだりします。


今回は以下のようになります。

被験者間要因：sex（A）
被験者内要因：W0,W1,W2の繰り返し（B）

今回のExcelではsexの列が1番右にありますが被験者間要因を左、被験者内要因を右に書くのが習わしのようなのでAsBデザインということになります。

被験者内要因（改善の傾向）に性別が影響しているかどうかを見たいということになります。

3.反復測定分散分析（repeated measures ANOVA）を行う

前回の記事で正規分布の確認を行いましたので、早速反復測定分散分析を行います。

反復測定分散分析を行うには主に2つの方法があります。

１．carパッケージのAnova関数を使う（EZRもcarパッケージを使ってる）
２．ANOVA君を使う

ANOVA君は井関龍太さんが作成した分散分析に特化したスクリプトです。
使い方に特徴がありますが使えるようになると慣れると簡単でどの方法よりも詳細に分析できます。

ここではまず1の方法を紹介します。

反復測定分散分析を行うには次の方法になります。

【Anova関数を使って反復測定分散分析】
①carパッケージを読み込む
②被検内要因（反復測定）は順番を指定するためのdata.frameを作る必要がありdata.farmeとfactor型のベクトルを事前に作っておく
③lm関数でモデルを作る（後のツイート参照）
④Anova関数で指定する（その後のツイート参照）

— MITTI1210 (@MITTI12101) November 19, 2019

①carパッケージを読み込む
ここは前回の記事と全く同じです。
まずはcarパッケージを読み込みます。
carパッケージはEZRをインストールする時に一緒にインストールされています。
もしインストールしていなければ先にインストールします。

install.packages("car")

その後library関数で呼び出します。

library(car)

②被験者要因を指定するdata.frameを準備する
ここは前回の内容と全く同じです。
今回のデータではsexの列は「MとFの2つのデータがある」とパソコンは読み取れますがW0,W1,W2はwideデータとなっているため「この3列が繰り返しのデータである」ことがパソコンは読み取れません。それがわかるためのものを2つ作ります。

まず繰り返しの順番がわかるtimeというベクトルを作ります。
このデータはfactor型である必要があるためfactor関数を使います。
そしてtimeのデータをdata.frame型に変更したidataを作ります。

 
そして大事なのがlevels = です。levels = を指定しないとfacorは五十音順になります。
今回は問題ないのですがもしpre,mid,postのように列名を付けるとlevels = を使わないとmid→post→preの順番になってしまいます。
factor関数に関しては【1-10】Rでよく使われる型について説明します。
 

③lm関数でモデルを作る

次に分析するためのモデルを作ります。ここが前回の記事と唯一違うところです。

モデルはlm関数を使います。model（名前は何でもいい）という変数名に結果を入れます。

被験者内要因は3列になっているのでcbind関数でつなげます。
被験者間要因は~の右に入れます。もし被験者間要因が複数あれば * でつなぎます（交互作用を見たい場合）。もし交互作用を見ない場合は + でつなぎます。

そして必ず大切なのがcontrasts = の箇所です。
数学的には反復測定分散分析はタイプⅢ平方和というものを計算します。
詳しい話は井関龍太さんの記事にありますが、contrasts = を設定しないと正しいタイプⅢ平方和を計算してくれないため、この設定が必要となります。
使い方はcontrasts = list(〇 = contr.sum, △ = contr.sum, □ = contr.sum)のようにlist内で指定します。

④Anova関数で分散分析を実行する

ここは前回の記事と全く同じです。
反復測定分散分析を行う時はcarパッケージのAnova関数を使います。Aは必ず大文字です。
res（resultの略、何でも良い）という変数名に結果を入れます。

Anova(lm関数で作ったモデル, idata = ◯◯, idesign = ~ ◯◯, type = "Ⅲ")

lm関数で作ったモデル：先程作ったmodel
idata：先程つくったidata
idesign：先程作ったtime（timeの前に ~ を付けます）
type："Ⅲ"または3とする（反復測定分散分析の場合は必ずⅢ）

⑤summary関数で結果を表示する

summary(Anova関数の結果, multivariate = FALSE)

multivariate = FALSEと入れるとEZRと同じ結果が表示されます。

①〜⑤をまとめると以下になります

library(car)
time <- factor(c("W0", "W1", "W2"), levels = c("W0","W1", "W2"))
idata <- as.data.frame(time)
model <- lm(cbind(W0,W1,W2) ~ sex, 
            data = grip_rep_anova, 
            contrasts = list(sex = contr.sum))
res <- Anova(model, idata = idata, idesign = ~time, type = 3)
summary(res, multivariate = FALSE)

行が長くて大変ですが実は色の部分以外はコピペ可能です。
（被験者間要因があればまた変わります。次回記事で紹介します）


結果はまずMauchlyの球面仮定を確認し、有意差がなければ上、あれば下をチェックします。

sex→握力は性別により変化する
time→握力は時間経過により変化する
sex:time→交互作用（時間経過に伴う握力の改善の程度は男女で異なるか）


結果はEZRと同じになりました。

4.まとめ
今回は2要因以上ある反復測定分散分析を行いました。
前回の記事との違いは被験者間要因があるかどうかでした。
被験者間要因があるとcontrasts = 設定が必要になりますので注意が必要です。

そして前回・今回はEZRと同じ方法のcarパッケージのAnova関数を使いました。
ただどうしてもコードが長くなってしまいます。
ただRにはのANOVA君という分散分析に特化した非常に有名なスクリプトがあります。

次回はANOVA君を使った分散分析を紹介します。

5.今回使ったコード

#データのダウンロード
url <- "https://haru-reha.com/wp-content/uploads/2018/05/demo-repeated-measures-anova.xlsx"
destfile <- "demo-repeated-measures-anova.xlsx" #今まで第4章の記事を読まれて実践されてる方は destfile <- "data/demo-repeated-measures-anova.xlsx"

download.file(url, destfile)

#データの読み込み
library(readxl)
grip_rep_anova <- read_excel("data/demo-repeated-measures-anova.xlsx", 
                             range = "B2:E32")
View(grip_rep_anova)

#sexの列もいれた反復測定分散分析
library(car)
time <- factor(c("W0", "W1", "W2"), levels = c("W0","W1", "W2"))
idata <- as.data.frame(time)
model <- lm(cbind(W0,W1,W2) ~ sex, 
            data = grip_rep_anova, 
            contrasts = list(sex = contr.sum))
res <- Anova(model, idata = idata, idesign = ~time, type = 3)
summary(res, multivariate = FALSE)