Programming
category
Programming: Rの覚書
| Rの覚書 |
- 多くのWebサイトを参考にさせていただきました。ごちゃまぜになっているので各サイト名を記すことはできませんが、みなさまありがとうございました。
- 自分の覚書きとして掲載しているため、正確でない記述が多くあります。
- 整理中のため、汚いまま残っています。
- 専門というわけではないので、良いコードではないものも多いことは理解しています。
| 目的 | コード |
| ライブラリ | library (ggplot2)#グラフ描画 library (grid) #unit関数を使うためにロード library(psych) #基本統計量の計算describe用 library(tidyr)# データの整形:ロング、ワイド library(dplyr)#group_by, library(devtools) #geom_pointのstroke設定を使うのに必要 library(lmerTest);#混合線形モデルlmeのため |
| 変数の全削除 | rm(list = ls(all.names = TRUE)) #変数の全削除 |
| データの読み込み: クリップボード | x <- read.table(“clipboard” ,header = T) |
| データの読み込み: txt形式、タブ区切り、ヘッダーTRUE | #データの読み込み df <- read.table(‘data.txt’ , skip = 0, header = T, sep = “\t”, stringsAsFactors = FALSE) |
| データの読み込み: csv形式 | #stringsAsFactors = FALSEにしないと、list型でグループ化されて、計算できない df <- read.table(‘data.csv’ , skip = 0, header = T, sep = “,”, stringsAsFactors = FALSE) |
| データセットの読み込み | df <- iris |
| ファイル名の取得、拡張子の削除 | #フォルダに含まれるtxtファイル名を取得 FileName <- dir(pattern = “\\.txt$”); #substrで拡張子を削除、ncharで文字数の取得 FileID <- substr(FileName, 1, nchar(FileName) – 4); |
| フォルダ内の特定種類のファイルを読み込んで、rbindで1つにまとめる | library(dplyr) #パイプ処理のため #csvファイル一覧を、lapplyでread.csvして、dfとして一括読み込み df <- dir(pattern = “\\.csv$”) %>% lapply(read.csv); #do.callでrbindを呼び出してdfに適用。rbindで1つのフレームにする。forより早い。 dfbind <- do.call(rbind, df);#列数のエラー確認用 for(i in 1:length(df)){ print(ncol(df[[i]]))} |
| フォルダ内の特定種類のファイルを読み込んで、rbindで1つにまとめる: skip やheaderの指定がある場合 | library(dplyr); #pipe処理用 #フォルダに含まれるcsvファイル名を取得 #csvファイル一覧のリストを、lapplyでread.csvして、dfとして一括読み込み #skipしないと、列数が3列になる。最初の5行で列数を判断しているため。 df <- dir(pattern = “\\.txt$”) %>% lapply(read.csv,sep=”,”,header=F, stringsAsFactors = FALSE,skip=27);#do.callでrbindを呼び出してHbDataに適用。rbindで1つのフレームにする。forより早い。 dfbind <- do.call(rbind, df);#列数のエラー確認用 for(i in 1:length(df)){ print(ncol(df[[i]]))} |
| エクセルexcelファイルの読み込み | library(readxl)#xlsの読み込み用パッケージ #col_typesで日付などのデータを書式設定 df <- dir(pattern = “\\.xls$”) %>% lapply(read_excel, col_types = c(“date”,”date”,”guess”,”guess”,”guess”,), skip=1); dfbind <- do.call(rbind, df);#列数のエラー確認用 for(i in 1: length(df)){ print(ncol(df[[i]]))} |
| フォルダ内の特定種類のファイル名を全取得して、処理を繰り返す | df1 <- dir(pattern = “\\.csv$”);#フォルダに含まれるTXTファイル名を取得 for (i in 1: length(df1)){#ファイル数だけ繰り返すfordf <- read.csv(df1[i], skip = 0, header = F, stringsAsFactors = FALSE) #stringsAsFactors = FALSEにしないと、list型でグループ化されて、計算できない}#ファイル数だけ繰り返すfor用 |
| 複数行のコメントアウト | #複数行のコメントアウト ココから—– if(0){ print(dfir) print(dfir) } #複数行のコメントアウト ココまで—– |
| 実行スクリプトのパスを取得し、作業ディレクトリにする | frame_files <- lapply(sys.frames(), function(x) x$ofile) frame_files <- Filter(Negate(is.null), frame_files) PATH <- dirname(frame_files[[length(frame_files)]]) |
| 文字の接続: スペース無し | paste(c(HbFile[i],”.png”), collapse=”” ) |
| 空フレームの作成 | Rate <- data.frame();#すべてのデータを格納する空フレーム |
| 変数の保存: csvファイル | write.csv(df,”Rate.csv”,quote=F, row.names=F) write.csv(df,paste(c(FileID[i],”.csv”), collapse=”” ),quote=F, row.names=F) |
| 連番の変数を作成する | for (i in 1:5) assign(paste(“A”, i, sep=””), i) # 変数 Ai に i を代入 |
| ベクトルの値を使って、変数を作成する ファイル名を変数名として指定する。 Assignで変数の作成 |
FileName <- dir(pattern = “\\.csv$”) #substrで拡張子を削除、ncharで文字数の取得 FileID <- substr( FileName, 1, nchar( FileName) – 4); FileNum <- length(FileNameu) for (i in 1:FileNum){ assign(FileID[i], read.table(FileName[i], skip = 1, header = F, stringsAsFactors = FALSE)) } |
| 複数の変数に、同様の処理をおこなう eval: Rスクリプトの実行 paste: 変数名を格納したベクトルと数式を接続 |
for (i in 1:FileNum){ eval(parse(text=paste(FileID[i],”[,4] <- “,FileID[i],”[,1]”,sep = “”))); } |
| 列名の変更: パイプ処理 | df <- esoph # esophデータセットの読み込み df2 <- df %>% `colnames<-` (c(“a”,”b”,”c”,”ch1″, “ch2”)) もしくは df2 <- df %>% setNames(.,c(“a”,”b”,”c”,”ch1″, “ch2”)) |
| データフレームへの列追加 | dfir <- iris[sample(1:nrow(iris),20),] #dfirにaddという列を追加する dfir2 <- transform(dfir,add=rep(1:5,4)); |
| 行・列の削除 | df <- iris df[,-5] #5列目を削除する |
| データフレームを縦につなげる（行の連結）: つなげるデータフレームの列名は、同じでなければならない。列の並び順が異なっていても、名前でおそらく結び付けられる。 | df1 <- ToothGrowth df2 <- ToothGrowth df3 <- rbind(df1,df2) |
| データフレームを横につなげる（列の連結）: つなげるデータフレームの列名は、同じでなければならない。列の並び順が異なっていても、名前でおそらく結び付けられる。 | df1 <- ToothGrowth df2 <- ToothGrowth df3 <- cbind(df1,df2) |
| データをロング型に整形する | library(tibble) #rownames_to_columnのため library(tidyr) #gatherのため library(dplyr) #パイプ処理のため mini_iris <- iris %>% group_by(Species) %>% slice(1) #1を1:2とかすると、1行目と2行目がスライスされる。Headとは異なる。 #複数列の値を1列にまとめたとき、まとめた列名をkey, その値列をvalueで指定する。マイナス指定することで、その列の値は表示されない。 mini_iris %>% gather(key = flower_att, value = measurement, -Species) #パイプ処理しない場合 mini_iris <- rownames_to_column(as.data.frame(mini_iris),’id’) gather(data = mini_iris, key = flower_att, value = measurement, -Species, -id) |
| データをワイド型に整形する | library(tibble) #rownames_to_columnのため library(tidyr) #gatherのため library(dplyr) #パイプ処理のため mini_iris <- iris %>% group_by(Species) %>% slice(1) # gatherで一度ロング型にする miniLong <- mini_iris %>% gather(key = flower_att, value = measurement, -Species) #spreadで、ワイド型に戻す。keyが、列名に、valueが各セルの値になる。 spread(data = miniLong,key = flower_att, value = measurement) |
| データフレームの転置 | df <- irisdf.t1 <- data.frame(t(iris)) #data.frameを宣言しないと、型がmatrixになる。 df.t2 <- t(iris)#matrixになることを確認class(df)class(df.t1)class(df.t2) |
| 条件 ifelse | dfir <- iris[sample(1:nrow(iris),20),] dfir[,6] <- ifelse(dfir[,2] > 3.4, 1, 0) |
| 値に特定文字を含むか否かの判定 | dfir <- iris[sample(1:nrow(iris),20),] ifelse(grepl(“g”, dfir$Species)==”TRUE”,1,0) |
| 条件のor (ベクトルの、対応する全要素を比較) | dfir <- iris[sample(1:nrow(iris),20),] dfir[,2] > 3.4 | dfir[,1] > 4.8 |
| 条件のor (ベクトルの最初の要素のみ比較) | dfir <- iris[sample(1:nrow(iris),20),] dfir[,2] > 3.4 || dfir[,1] > 4.8 |
| 条件のand (ベクトルの、対応する全要素を比較) | dfir <- iris[sample(1:nrow(iris),20),] dfir[,2] > 3.4 & dfir[,1] > 4.8 |
| 条件のand r (ベクトルの最初の要素のみ比較) | dfir <- iris[sample(1:nrow(iris),20),] dfir[,2] > 3.4 && dfir[,1] > 4.8 |
| 条件で行の抽出 | library(dplyr) #カンマで連結すれば&条件 dfir <- iris[sample(1:nrow(iris),20),] %>%filter(Species == “setosa” , Petal.Width == 0.2)#パイプ処理しない場合 dfir <- iris[sample(1:nrow(iris),20),]filter(dfir ,Species == “setosa” , Petal.Width == 0.2)#|で連結すればOR条件 dfir <- iris[sample(1:nrow(iris),20),]filter(dfir ,Species == “setosa” | Petal.Width == 1.8) |
| 条件で行の抽出 | library(dplyr) dfir <- iris[sample(1:nrow(iris),20),] %>%subset( grepl(“se”, Species)); |
| 等しくない値を抽出。 | # != が～でないを示す。 dfir <- iris[sample(1:nrow(iris),20),] dfir$Species != “setosa” |
| 欠損値でない値を抽出 | dfir <- iris[sample(1:nrow(iris),20),] dfir[sample(1:20,7),3] <- NA; # !is.naが、NAでない、を示す。 dfir[!is.na(dfir[,3]),] |
| 無作為にサンプリングする。sample関数 | #無作為サンプリング。1-150の中から100個を、重複なしにサンプリングする。 sample(1:150,100,F) #1-10の中から、100個を重複ありでサンプリングする。 sample(1:10,100, T) #1-10の中から、100個を重複なしではサンプリングできないので、エラーが出る。 sample(1:10,100, F) |
| 特定の内容を繰り返す | rep(1:5,4) |
| ある範囲を公差で分割する | seq(0,8,by =0.6)# seq(a, b, by = c)a から b まで c づつ増加するベクトルを生成． |
| チャネル数などの名称を繰り返すベクトル作成 | c(“Cond”, “Subject”, “Trial”,”Time”, paste(“ch”,rep(seq(1,4,1), each=3),rep(c(“oxy”,”deoxy”,”total”), 4),sep = “”)) |
| 列ごとに計算する | df <- iris #第2項が2.文字が含まれる列が1列でもあるとエラーが出る。 apply(df[,-5],2,mean) #1回目Restの平均値 |
| 行ごとに計算する | df <- iris #第2項が1.文字が含まれる行があるとエラーが出る。 apply(df[,-5],1,mean) #1回目Restの平均値 |
| 基本統計量の計算：標準誤差SEを含む | library(psych) # describeByのためのパッケージ df <- ToothGrowth describeBy(df) |
| グループ化して集計(1) 複数列について、集計の値をもとめる |
df <- iris #複数列について、1列の値を基にグループ化して集計 #グループ化する列は、listに要変換 aggregate(x=df[c(“Sepal.Length”, “Sepal.Width”)], by=list(species=df$Species), FUN=mean) |
| グループ化して集計(2) 複数列でグループ化して、集計する |
df <- ToothGrowth # ToothGrowthデータセットの読み込み aggregate(x=df[“len”], by=list(supp=df$supp, dose = df$dose ), FUN=mean) #ncasesとncontrolsを、alcgpとtobgpでグループ化して、meanとsumを出す。 aggregate(cbind(ncases, ncontrols) ~ alcgp + tobgp, data = esoph, FUN = function(x) c(mn = mean(x), n = sum(x) )) |
| グループ化して集計(3) 複数列でグループ化して、集計する |
#こっちがオススメ。SEも出る。 library(psych) # describeByのためのパッケージdescribeBy(sat.act$age,list(sat.act$education,sat.act$gender), mat=TRUE)df <- Moore # Mooreデータセットの読み込み # conformityとfscoreを、partner.statusとfcategoryでグループ化して、meanとsumを出す。aggregate(cbind(conformity, fscore) ~ partner.status + fcategory, data = df, FUN = function(x) c(mn = mean(x), n = sum(x) )) |
| グループ化して集計(4) 複数列でグループ化して、関数を複数使って集計する |
library(dplyr) #パイプ処理と、summarize関数のためのlibrary df <- ToothGrowth #suppとdoseでグループ化して、lenのmeanとsumを求める df %>% group_by(supp, dose) %>% summarize(mean = mean(len),sum =sum(len)) |
| 2変量のクロス集計 | df <- ToothGrowth table(df[,2:3]) |
| 3変量のクロス集計 | df <- ToothGrowth #c関数内の並び方で、集計方法が変わる。Lenとsuppでクロス集計している df2 <- with(df, table(df[,c(“dose”,”len”,”supp”)])) ftable(df2) |
| 等分散を仮定しない Welch’S t-test | #等分散を仮定しない Welch’S t-test df <- PlantGrowtht.test(df[df$group==”ctrl”,1], df[df$group ==”trt1″,1], var.equal=F) |
| matlabのrepmat関数を導入 | #matlabのrepmat関数を導入——————– #http://www.waxworksmath.com/Authors/G_M/Hastie/Code/Chapter4/repmat.R repmat = function(X,m,n){ ##R equivalent of repmat (matlab) mx = dim(X)[1] nx = dim(X)[2] matrix(t(matrix(X,mx,nx*n)),mx*m,nx*n,byrow=T) } |
| plot用 | |
| 描画キャンパスのレイアウト | par(mfrow=c(3,2))# 描画キャンパスを6分割する par(mfrow=c(1,1))# 描画キャンパスを1分割に戻す |
| plotでの描画 | dfir <- iris #Ch1折れ線グラフの作成。C: X、Y座標を指定。ylim:Y軸の幅、colで色、lwd: 太さ plot(dfir[,c(3,2)],ylim = c(0,10),type=”l”,col=1, lwd =5 ) lines(dfir[,c(3,1)],type=”l”,col=2, lwd =1) |
| グラフを重ねて描く: 線グラフ Plotしたあとに、linesで上描きできる。 |
dfir <- iris plot(dfir[,c(3,2)],ylim = c(0,10),type=”l”,col=1, lwd =5 ) lines(dfir[,c(3,1)],type=”l”,col=2, lwd =1) |
| グラフを重ねて描く: 散布図 Plotしたあとに、pointsで上描きできる。 |
dfir <- iris plot(dfir[,c(3,2)],ylim = c(0,10),type=”l”,col=1, lwd =5 ) points(dfir[,c(3,1)] , pch =21, col=”#0070c0″, bg =”#ffffff”, lwd =2, cex = 2) #pch:プロットする図形の形 #cex: プロットサイズ #bg: fillの色。ただし、pch=21-25にしか適用できない。 #PointShape = c(21,16,24,4,16,16) #散布図の形の順番指定。〇、●、△、×、● #色の設定（line用）（青、ｵﾚﾝｼﾞ、黒、,明るいグレー,グレー） #colours_line <- c(“#0070c0”, “#e46c0a”, “#000000″,”#c8c8cb”, “#7f878f”) #色の設定（fill用）（白、ｵﾚﾝｼﾞ、黒、,明るいグレー,グレー） colours_fill1 <- c(“#ffffff”, “#e46c0a”, “#000000”, “#c8c8cb”, “#7f878f”) |
| グラフを重ねて描く: par(new=T)#上書き指定 X・Y軸座標値、ラベルなどの指定が必要なので、要注意。 |
dfir <- iris plot(dfir[,c(3,2)],ylim = c(0,10), xlim = c(0, 10), type=”l”,col=1) par(new=T)#上書き指定 #xlabとylabをブランクにする。また、ylimとxlimも同じに設定しないと、文字がズレる。2枚のグラフを重ね描きしているので、それぞれの文字も2回描かれているため。 plot (dfir[,c(3,1)] ,ylim = c(0,10), xlim = c(0, 10), type=”both”,col=2, ylab=””, xlab=””) |
| ggplot2用 | |
| 色の設定（line用）（青、ｵﾚﾝｼﾞ、黒、,明るいグレー,グレー） | #色の設定（line用）（青、ｵﾚﾝｼﾞ、黒、,明るいグレー,グレー） colours_line <- c(“#0070c0”, “#e46c0a”, “#000000″,”#c8c8cb”, “#7f878f”) |
| 色の設定（fill用）（白、ｵﾚﾝｼﾞ、黒、,明るいグレー,グレー） | #色の設定（fill用）（白、ｵﾚﾝｼﾞ、黒、,明るいグレー,グレー） colours_fill1 <- c(“#ffffff”, “#e46c0a”, “#000000”, “#c8c8cb”, “#7f878f”) |
| 散布図の図形 | PointShape = c(21,16,24,4,16,16) #散布図の形の順番指定。〇、●、△、×、● |
| 使えるフォントの確認 | names(windowsFonts()) #使えるフォントの確認 |
| 使うフォントの設定: メイリオ | windowsFonts(MEI = windowsFont(“Meiryo”)) |
| 使うフォントの設定: MigMix 1P | windowsFonts(Mig = windowsFont(“MigMix 1P”)) |
| グラフを画像として保存 #R上にいちいち描画せず、直接出力するので、保存が早いはず | png(paste(c(HbFile[i],”.png”), collapse=”” ), width = 800, height = 450,pointsize = 30) # 描画デバイスを開く par(lwd = 2,bg =”black”)#線の太さと背景色 plot(g)dev.off() |
| jpeg(paste(“temp”,”.jpg”), collapse=”” ), width = 1280, height = 720,pointsize = 30, quality = 100) ;#出力するときの文字サイズ。ここで設定する必要あり。 par(lwd = 2,bg =”black”)#線の太さと背景色 |
|
| Z-scoreの計算 | # Z-scoreでの分析開始————————– #文献：2015NIRSを用いた局所温冷刺激に対する脳賦活反応解析 Zbase <- data.frame(PreMean = apply(HbData[ZscoreS:ZscoreL,5:24], 2, mean),#Pretaskの平均 PreSd = apply(HbData[ZscoreS:ZscoreL,5:24], 2, sd));#PretaskのSD #補正後のデータ用変数HbDataR HbDataZscore <- HbData; #標準得点(Z-score)化 for (m in 1:20){ HbDataZscore[,4+m] <- (HbDataZscore[,4+m]-Zbase[m,1])/Zbase[m,2]#補正 } |
| 移動平均 | #移動平均のための時系列分析——————————- #時系列オブジェクトの作成 IIP <- ts(HbDataR, start=0, end =120, frequency=1/0.015); #ts関数が時系列オブジェクト #http://www.is.titech.ac.jp/~mase/mase/html.jp/temp/ts.jp.html#filter関数により移動平均を計算 w5 <- rep(1,334)/334;#5秒平均#filterは移動平均の各項の重みを与える引数 #methodは”convolution”とすれば移動平均、”recursive”とすれば指数平滑化 #filter=w3の場合、side=2だと(yt-1+yt+y+1)/3、side=1だと(yt+yt-1+yt-2)/3となる#IIP1 <- filter(IIP,filter=w1,method=”convolution”,sides=2) IIP5 <- filter(IIP,filter=w5,method=”convolution”,sides=2);#matlabのrepmat関数を導入——————– #http://www.waxworksmath.com/Authors/G_M/Hastie/Code/Chapter4/repmat.R repmat = function(X,m,n){ ##R equivalent of repmat (matlab) mx = dim(X)[1] nx = dim(X)[2] matrix(t(matrix(X,mx,nx*n)),mx*m,nx*n,byrow=T) } #60sで0になるように補正—————————-#60sで0になるように補正する変数IIP5r IIP5r <- IIP5; #write.csv(IIP5,”IIP5.csv”,quote=F, row.names=F) #write.csv(IIP5r,”IIP5r.csv”,quote=F, row.names=F)#60sで0になるように補正 IIP5r[,5:24] <- IIP5r[,5:24]-repmat(matrix(IIP5r[4001,5:24],1,20),nrow(IIP5r),1)#—————— IIP5r[,1] <- Cond[i];#条件 IIP5r[,2] <- Subject[i];#被験者ID IIP5r[,3] <- Trial[i];#Trial回数#補正したデータをcsvで書き出し————————–#assignでFileIDを変数名とし、補正を終了したHbDataを格納する #assign(FileID[i], HbData)write.csv(IIP5r[4001:6001,],paste(c(FileID[i],”.csv”), collapse=”” ),quote=F, row.names=F)#Z-scoreでの分析開始————————– #文献：2015NIRSを用いた局所温冷刺激に対する脳賦活反応解析 Zbase <- data.frame(PreMean = apply(HbData[ZscoreS:ZscoreL,5:24], 2, mean),#Pretaskの平均 PreSd = apply(HbData[ZscoreS:ZscoreL,5:24], 2, sd));#PretaskのSD#補正後のデータ用変数HbDataR HbDataZscore <- HbData;#標準得点(Z-score)化 for (m in 1:20){HbDataZscore[,4+m] <- (HbDataZscore[,4+m]-Zbase[m,1])/Zbase[m,2]#補正 }#移動平均のための時系列分析——————————- #時系列オブジェクトの作成 IIPZscore <- ts( HbDataZscore, start=0, end =120, frequency=1/0.015); w5 <- rep(1,334)/334;#5秒平均 IIP5Zscore <- filter(IIPZscore,filter=w5,method=”convolution”,sides=2);#60sで0になるように補正する変数IIP5r IIP5rZscore <- IIP5Zscore; #60sで0になるように補正 IIP5rZscore[,5:24] <- IIP5rZscore[,5:24]-repmat(matrix(IIP5rZscore[4001,5:24],1,20),nrow(IIP5r),1) |
| ##inputデータは経済産業省より取得 data <- read.csv(“http://www.meti.go.jp/statistics/tyo/iip/result/h2afdldj/csv/ha2nsgo3j.csv”,header=TRUE)#データの整形 delete <- c(1,3) data <- data[-delete,] colnames(data) <- data[1,]new_name <- NULL for (i in 1:ncol(data)) { new_name <- c(new_name,as.character(data[1,i])) }names(data) <- new_name data <- data[-1,] names(data)[1] <- “年月” head(data)#時系列データに変換 data <- data[1:359,] IIP <- ts(data, start=c(1978,1),frequency=12) #ts関数が時系列オブジェクトhttp://www.is.titech.ac.jp/~mase/mase/html.jp/temp/ts.jp.html ts.plot(IIP[,2],col=c(1:3))#filter関数により移動平均を計算 w3 <- rep(1,3)/3; w7 <- rep(1,7)/7; w13 <- rep(1,13)/13 w25 <- rep(1,25)/25; w49 <- rep(1,49)/49 w98 <- rep(1,98)/98w3 <- rep(1,3); w7 <- rep(1,7); w13 <- rep(1,13) w25 <- rep(1,25); w49 <- rep(1,49) w98 <- rep(1,98)IIP3 <- filter(IIP,filter=w3,method=”convolution”,sides=2) IIP7 <- filter(IIP,filter=w7,method=”convolution”,sides=2) IIP13 <- filter(IIP,filter=w13,method=”convolution”,sides=2) IIP25 <- filter(IIP,filter=w25,method=”convolution”,sides=2) IIP49 <- filter(IIP,filter=w49,method=”convolution”,sides=2) IIP98 <- filter(IIP,filter=w98,method=”convolution”,sides=2) #filterは移動平均の各項の重みを与える引数 #methodは”convolution”とすれば移動平均、”recursive”とすれば指数平滑化 #filter=w3の場合、side=2だと(yt-1+yt+y+1)/3、side=1だと(yt+yt-1+yt-2)/3となる par(mfrow=c(3,2)) #図を6分割する ts.plot(IIP[,5],type=”l”,main=”原系列”) ts.plot(IIP3[,5],type=”l”,main=”3期移動平均”) ts.plot(IIP13[,5],type=”l”,main=”13期移動平均”) ts.plot(IIP25[,5],type=”l”,main=”25期移動平均”) ts.plot(IIP49[,5],type=”l”,main=”49期移動平均”) ts.plot(IIP98[,5],type=”l”,main=”98期移動平均”)par(mfrow=c(1,1)) #図を1分割に戻す |
|
| 分散分析 | source(“anovakun_481.txt”) #分散分析をするときには、factorにすることを忘れない AkunDfNS <- data.frame(s = factor(TotalAno$Sub), a = factor( TotalAno$Ex), b = factor(TotalAno$PrePost), result = TotalAno$nau)#2要因とも対応がある2要因分散分析(RBFpqデザイン) (p.116) #print(summary(aov(result ~ a * b + Error(s + s:a + s:b + s:a:b), AkunDfNS)))#getaが効果量詳細は、ANOVA君/反復測定デザインにおける効果量 #cindが信頼区間,cindは平均値の差、cipairは多重比較に興味があるとき #詳細はANOVA君/反復測定デザインにおける信頼区間 anovakun(AkunDfNS[1:56,], “sAB”, 2,2, long = T,geta = T,cind = T,cipair = T) a <- summary(anovakun(dat, “sABC”, 2,16,6, long = T,geta = T,cind = T,cipair = T)) write.csv(anovakun(dat, “sABC”, 2,16,6, long = T,geta = T,cind = T,cipair = T),”FacScore.csv”,quote=F, row.names=F) |
| 混合線形モデルの計算 | #混合線形モデルの計算 m1 <- lmer(result~a*b + (1|s) + (1|s:a) + (1|s:b), data = AkunDfNS[1:56,]); anova(m1, ddf = “Kenward-Roger”); jigo <- step(m1, ddf = “Kenward-Roger”);#交互作用後の事後検定 |
| 相関係数、検定力、相関係数をz変換 | library(pwr); #相関係数の効果量の計算用 x <- c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1); y <- c( 2.6, 3.1, 2.5, 5.0, 3.6, 4.0, 5.2, 2.8, 3.8); #相関があるかどうかの検定pearsonの方法。他に2つある。 #”kendall” :ケンドールの順位相関係数の無相関検定を行う． #”spearman” ： スピアマンの順位相関係数の無相関検定を行う． result <- cor.test(x, y, method = “pearson”) #検定力の計算 #r相関係数、sig.level有意水準、n標本数 pwr.r.test(r=result$estimate,n=length(x),sig.level=result$p.value,alternative=”two.sided”)$power #z変換 1/2*log((1+ result$estimate)/(1- result$estimate)) |
| 混合線形モデルの計算 | #混合線形モデルの計算 library(lmerTest);#混合線形モデルlmeのため YoiData <-data.frame(s = factor(YoiData$Sub), a = factor(YoiData$Ex), b = factor(YoiData$YoiHIY), YoiData[,c(3:6,12:27)]); YoiDataR <-data.frame(YoiData[,1:4])m1n <- lmer(nau~a*b + (1|s) + (1|s:a) + (1|s:b), data = YoiData, control=lmerControl(check.nobs.vs.nlev = “ignore”, check.nobs.vs.rankZ = “ignore”, check.nobs.vs.nRE=”ignore”)); anova(m1n, ddf = “Kenward-Roger”); |
| データを集める | AveTime <- data.frame(matrix(nrow =1,ncol = 1)); for (i in 1:FileNum){#ファイル数だけ繰り返すfor#TXTファイルの読み込み.ファイル名を変える #skipで、上から35行を読み込まないと指定 #列名は設定しない #stringsAsFactors = FALSEにしないと、list型でグループ化されて、計算できない HbData <- read.table(HbFile[i], skip = 0, header = F, stringsAsFactors = FALSE, sep =”,”); s <- paste(“AveTime <- transform(AveTime,”, Subject[i], ” = HbData[1:74810,2])”, sep=””) eval(parse(text=s)) #AveTime <- transform(AveTime,a1 = HbData[1:74810,2] )}#ファイル数だけ繰り返すfor用 |
| 線形混合モデル | library(lmerTest);#線形混合モデルlmerのため library(MuMIn);#線形混合モデルのR2の計算のため。r.squaredGLMM関数#混合線形モデルの計算 YoiData <-data.frame(s = factor(YoiData$Sub), a = factor(YoiData$Ex), YoiData[,3:18]);#——-Ch1———————- m1z1in <- lmer(z1~a + (1|s) , data = YoiData); summary(m1z1in) r.squaredGLMM(m1z1in)#GLMM用のr2 anova(m1z1in, ddf = “Kenward-Roger”); |
