色々と分析する中であんなことやりたい、こんなことやりたい、というのが出てくるのだが、覚えきれずに毎回調べる羽目になるので、ここで頑張ってつくってみる。今後、少しずつ増やしていく。
相関行列
標準パッケージのstats::corかpsychパッケージのpsych::corr.testか。
corだと、NAを含んだデータの場合、useオプションを適切に設定しないとNAが返されるので注意。useオプションの設定値はRのhelp(cor)参照。
stats::cor(data, use="pairwise.complete.obs")
psych::corr.test(data)
相関行列を図示することもできる。詳しくは以下のリンク。
NAを含むデータの平均
NAを含んでいるデータをそのまま平均だそうとするとNAしか返ってこない。
そういう時はmean関数のna.rmをTUREにする。
data<-c(1,2,3,NA,4,5,6)
mean(data) #NAが返る
mean(data, na.rm=TRUE)
グループ単位での平均値の求め、列に追加する
マルチレベル分析でよく使う。グループ単位でネストされたデータがあるとき、グループ単位での平均値を各行に追加したい。そういう時はaveを使う。
df <-data.frame(
grpid =c(1,1,1,1,1,2,2,2,2,2,3,3,3,4,4,4,4),
data =c(1,2,3,4,5,3,2,3,4,4,2,1,3,5,5,7,9)
)
ave(df$data, df$grpid)
# 出力：3.0 3.0 3.0 3.0 3.0 3.2 3.2 3.2 3.2 3.2 2.0 2.0 2.0 6.5 6.5 6.5 6.5
# 元のデータに追加する場合は以下。
df <-cbind(df,平均=ave(df$data, df$grpid))
もしデータにNAが含まれていた場合、そのままやると、そのグループの平均値はNAになる。
そういう場合はこう。
df <-data.frame(
grpid =c(1,1,1,1,1,2,2,2,2,2,3,3,3,4,4,4,4),
data =c(1,2,3,4,5,3,2,NA,4,4,2,1,3,5,5,7,9)
)
ave(df$data, df$grpid)
# 出力： 3.0 3.0 3.0 3.0 3.0 NA NA NA NA NA 2.0 2.0 2.0 6.5 6.5 6.5 6.5
ave(df$data, df$grpid, FUN=function(x) mean(x, na.rm=TRUE))
# 出力： 3.0 3.0 3.0 3.0 3.0 3.25 3.25 3.25 3.25 3.25 2.0 2.0 2.0 6.5 6.5 6.5 6.5
FUNでいろいろと関数を変えれて、ここにsd()やらmedian()なんかも書き込むことができる。
グループ単位での平均の一覧取得
psychパッケージのdescribeByでもできるが、以下のようにもできる。describeByだと不要なデータも出てくるし、グループごとのデータの取り出しも実は結構面倒（参考ページ）。なので、データ処理のために使うなら、以下の方が使い勝手が良い。
df <-data.frame(
grpid =c(1,1,1,1,1,2,2,2,2,2,3,3,3,4,4,4,4),
data =c(1,2,3,4,5,3,2,NA,4,4,2,1,3,5,5,7,9)
)
aggregate(data~grpid,df, FUN=mean)
# grpid data
# 1 1 3.00
# 2 2 3.25
# 3 3 2.00
# 4 4 6.50
aggregateにはほかにも書式があるので、Rのhelp機能からstats::aggregateを確認すべし。
Excelのvlookupみたいなことをする
各行のgrpidに対応する値を、別のデータフレームから検索して、対応する値を取得して、新しい列に入れこむというようなことをしたい。（例えば、grpidごとの部下データの平均値を、同じgrpidが振られた上司データに入れていく、みたいな）
そういう場合は、grepを使ってこうする。
(参考:https://sites.google.com/site/shichijostudents/home/r)
df.boss <-data.frame(
grpid =c(2,3,1,4,4,2,3,1,5)
)
df.subgroup <-data.frame(
grpid =c(1,2,3,4),
group_m = c(3.2,2.1,4.0,3.3)
)
group_m <- rep(NA, nrow(df.boss))
for( i in 1:nrow(df.subgroup) ){
target <- grep(df.subgroup$grpid[i], df.boss$grpid)
group_m[target] <- df.subgroup $group_m[i]
}
df.boss<-cbind(df.boss, group_m)
回帰分析系のパラメータの信頼区間
lmやglmあるいはlmerなどモデル分析を行う場合の出力の信頼区間を求めたい場合にはconfintを使う。
x <-c(1,1,1,1,1,2,2,2,2,2,3,3,3,4,4,4,4)
y <-c(1,2,3,4,5,3,2,NA,4,4,2,1,3,5,5,7,9)
df <-data.frame(x,y)
res<-lm(y~x,df)
# summary(res)
# ~~~~略~~~~
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 1.5115 1.0687 1.414 0.179
# x 0.9425 0.4039 2.333 0.035 *
# ---
confint(res)
# 2.5 % 97.5 %
# (Intercept) -0.78054596 3.803534
# x 0.07621897 1.808839
confintにはほかにもパラメータを設定できる。
回帰分析系の結果の図示と表示
sjPlotパッケージを使うのが良き。
lmでも,lmerでも、glmでも対応してくれる。
library(sjPlot)
res <- lmer( y~x1*x2+x3+(1|group), data)
plot_model(res,"pred",term=c("x1","x2") #交互作用項の単純傾斜分析
tab_model(res)
また日本語の対応もうまく行ってない。
さらにHTML出力やtexファイルの出力は出来るが、PDFへのknitには対応していない様子。
PDFへのknit出力するには、いくつか方法が提案されている。
ただ、どの方法もうまく行かなかった・・・。
探索的因子分析
スクリープロット、MAP基準やBIC基準など、んで、FA。FAの数字は抽出する因子数。回転については、デフォルトはオブリミン回転。Promaxを使おうと思うとGPArotationパッケージが必要。
library(psych)
library(GPArotation)
fa.parallel(data, fa = "fa", use = "pairwise.complete.obs")
vss(data)
fa(data,5,fm="ml",rotate="promax")
Lavaanを使った構造方程式モデリング関係
モデルの表記やら適合度指標の取得やら標準化解やら。
以下にまとめてる。
マルチレベルSEM
モデル式をlevel:で分けるのと、sem()にグループ識別子をcluster=で与えてやる。
ML.model <-'
Level:1
f1 =~ x1 + x2
f2 =~ x3 + x4
f1 ~ f2
Level:2
f1 =~ x1 + x2
f2 =~ x3 + x4
f1 ~ f2
'
fit <- sem(ML.model, data, cluster="grpid")
出力結果をいい感じの表にする
semoutputというパッケージを使う。CRANサーバからはインストールできない。インストールは以下。使う際には普通にlibraryで読み込む。
devtools::install_github("dr-JT/semoutput")
# ・・・
# lavaan でsemをする。結果はfitに保存。
sem_tables(fit)
なお、Rmarkdownで出力したときにHTMLのコードがそのまま表示されてしまう場合には、chankoptionに results='asis'を書き込む。
chankoptionの参考：
https://teramonagi.hatenablog.com/entry/20130615/1371303616
ただし、現時点ではMultilevel SEMの結果に対しては対応していない模様。
パス図を自動的に作る
semPlotパッケージとsemptoolsパッケージを使う。
レイアウトを変える
レイアウトをグリッドとして考えて、書き込みたい変数の位置を示すマトリクスを作る
library(semPlot)
library(semptools)
#・・・ x1 x2 x3 という外生変数とf1という内生変数があるモデルとする
# f1 を上に、x1～x3を下に書きたい.グリッドのc(行,列)でそれぞれの位置を指定
m <- layout_matrix(
f1 = c(1,2),
x1 = c(2,1),
x2 = c(2,2),
x3 = c(2,3)
)
# あるいは以下のようにしてもよい
m1<- matrix(c(c(0,"f1",0),c("x1","x2","x3")),3,2,byrow=F) # 列ごとに書く場合。要するにc(1列目の配置,2列目の配置,...)
m2<- matrix(c(c(0,"x1"),c("f1","x2"),c(0,"x3")),3,2,byrow=T)#行ごとに書く場合。要するにc(1行目の配置,2行目の配置,...)
semPaths(fit, "est", layout = m)
注意点として、Matrixに書き込む変数名はLavaanのモデルに書き込んだ変数名ではなく、semPathsが出力するノードのラベルを書かないといけない。つまり、モデル内で長い変数名をつけていても、semPathsのデフォルトでは3文字分しか出力されない。そしてMatrixも3文字に制限された変数名を書かないといけないということ。
長い変数名をそのまま書き込みたい場合には、semPathsの引数にnCharNodes=0を与えておく。そうするとモデルでつけてる変数名をMatrixにそのまま含めることができる（ただし、あまりに長い変素名をつけるとsemPathsが描くノードオブジェクトの変数名が非常に小さな文字になってしまって見にくい。
有意かどうかのアスタリスクをつける
SEMの結果から、推定値と有意確率をもとに、図の矢印に書き込むラベルを予め作成する。
# lavaan でのSEMの結果はfitに保存。
# 矢印に書き込むラベルを作成。ここでアスタリスクをつける。
table <- parameterEstimates(fit)
labels <-rep(NA,nrow(table))
for(i in 1:nrow(table)){
if(table$pvalue[i]<0.001) p[i]<-sprintf("%.2f,***",table$est[i])
else if(table$pvalue[i]<0.01) p[i]<-sprintf("%.2f,**",table$est[i])
else if(table$pvalue[i]<0.05) p[i]<-sprintf("%.2f,*",table$est[i])
else if(table$pvalue[i]<0.1) p[i]<-sprintf("%.2f,†",table$est[i])
else p[i]<-sprintf("%.2f",table$est[i]) #非有意なパスには係数を書きたくない場合にはこの行を消す
}
semPaths(fit, "est",edgeLabels = labels)
semptoolsを使って以下のようにも書ける。こちらの方がシンプル。alphaで基準と書き込むシンボルを設定。例えばn.s.を書き込みたくない場合には、"(n.s.)" = 1.00,を消せばよい。
ただし、例えば非有意のパスには推定値を書き込みたくない、といった細かな設定をしたい場合には、上の例の方法になる。
semPaths_plot <- semPaths(fit, "est",edgeLabels = labels, DoNotPlot=TRUE)
plot(p_pa<-mark_sig(semPaths_plot,fit,
alpha = c("(n.s.)" = 1.00, "†" = 0.10, "*" = .05, "**" = .01, "***" = .001)
)
)
なお、上の例でsemPaths()のなかのDoNotPlotは、semPathsは実行するけど、結果はPlotしないという命令。これをしておかないと、そのあとのmark_sigとで2つの図が表示されてしまう。
RmarkdownでPDF作成
こちらにまとめている。
Rmarkdownでタブ表記
## ○○○ の後に {.tabset}とつけると、そのセクションに属する下位セクションはタブ表記になる。