Hatena::ブログ(Diary)

arupaka-_-arupakaの日記

2017-01-20

Rで行列の文字コード変換

rownames(r_d)<-iconv(rownames(r_d),to="utf-8",from="cp932")
colnames(r_d)<-iconv(colnames(r_d),to="utf-8",from="cp932")

R3.00しかないオフライン環境でのlavaan(分散共分散解析 sem)のインストール

Step1: ネットのつながるパソコンで

 wget https://cran.r-project.org/src/contrib/Archive/lavaan/lavaan_0.5-17.tar.gz

0.5-17がR 3.0で使える一番新しいもの

Step2: インストールしたいパソコンにコピーして以下を実行

 R CMD INSTALL *0.5*17*tar.gz

ERROR: dependencies ‘mnormt’, ‘pbivnorm’, ‘quadprog’ are not available for package ‘lavaan

とでる。

Step3:

そこでネットのつながるパソコンでライブラリをゲット

wget https://cran.r-project.org/src/contrib/mnormt_1.5-5.tar.gz

wget https://cran.r-project.org/src/contrib/pbivnorm_0.6.0.tar.gz

wget https://cran.r-project.org/src/contrib/quadprog_1.5-5.tar.gz

Step4: これらをインストールしたいPCにアップロード

user1@panasonic ~

scp mnormt*tar.gz xxx@xxx:

scp pbivnorm*tar.gz xxx@xxx:

scp quadprog*tar.gz xxx@xxx:

Step5: インストール先でパッケージをインストール

R CMD INSTALL mnormt*tar.gz

R CMD INSTALL pbivnorm*tar.gz

R CMD INSTALL quadprog*tar.gz

Step6 lavaanをインストール

R CMD INSTALL lavaan*0.5*17*tar.gz

例:

library(lavaan)

x<-rnorm(100)

y<-3*x+rnorm(100,0,0.1)

z<-cbind(x,y)

model<-'y~x'

fit<-sem(model=model,data=z,estimator="ML")

summary(fit)

出力

lavaan (0.5-17) converged normally after 1 iterations

Number of observations 100

Estimator ML

Minimum Function Test Statistic 0.000

Degrees of freedom 0

Minimum Function Value 0.0000000000001

Parameter estimates:

Information Expected

Standard Errors Standard

Estimate Std.err Z-value P(>|z|)

Regressions:

y ~

x 3.009 0.010 301.411 0.000

Variances:

y 0.008 0.001

2017-01-18

Rで因子分析と分散共分散分析

Rで因子分析と分散分散分析をする。

  • 因子分析は、潜在要因をみつける方法。
    • 主成分分析と回帰分析と異なり観測変数xを説明する潜在変数をみつける方法(下のx1とx2を見つける)

y1=a*x1+a*x2

y2=a*x1+a*x2

y3=a*x1+a*x2

y4=a*x1+a*x2

    • 心理学等では、「瞬発筋力x1」や「持久力x2」や「器用さx3」など観測できない構成概念を扱うことが重要なためこの方法が発達した。
      • 観測できるのは、持久走のタイムy1、短距離走のタイムy2、ソフトボールなど記録y3、体前屈の記録y4、身長y5、体重y6など潜在要因からきまってくるもの
    • ただし、係数の決め方や変数の決め方に不定性(同じパラメータでおなじ説明力)があるので、軸の回転をおこなうことで、説明しやすい解や理解しやすい解をさぐる。
      • 各潜在要因が直交しない解がだすことができるのも因子分析の特徴
      • 自由に色々選べるので調整がききやすい
      • 代表的な自動回転法 promax回転(非直行)やvarimax回転(直行)。構造を目指す(出来る限り観測変数が一つの説明の要因で説明できるようにする回転。潜在要因と観測要因の関係が明確になるので説明しやすい)。
      • 要因の説明力をたかめると、個別データをバラバラにする野の違いもある(主成分と因子分析)

参考:Rで因子分析 商用ソフトで実行できない因子分析のあれこれ

http://www.slideshare.net/simizu706/r-42283141

#最低値がでる。BICがおすすめらしい

map <- vss(rr2, fm="minres")

#最大値がでる。BICがおすすめらしい

para <- fa.parallel(rr2)

        • ほかには主成分析→おりまがり+-1を調べる plot(eigen(cor(X))$values)
      • Step3 因子分析

#例:8要因でoblimin回転(promaxの進化版?)

v<-fa(rr3,nfactor=8,rotate = "oblimin")

#例:8要因でvarimax回転(直行)

v2<-fa(rr3,nfactor=8,rotate = "varimax")

#例:8要因でpromax回転(非直行)

v3<-fa(rr3,nfactor=8,rotate = "promax")

      • Step4: 結果の取り出し(各潜在要因のうち因子スコアが大きい観測量をとりだし)。

for(i in 1:8){

cat(labels*12[i],"\n")

s<-(v2$loadings[,i][rev(order(abs(v2$loadings[,i])))][1:10])

print(s)

#print(paste(names(s),collapse=", "))

#print(paste(sprintf("%.3f",s),collapse=", "))

}

#install.pakcages("lavaan")

#install.packages("semPlot")

#install.packages("knitr")

#install.packages("data.table")

#library(semPlot)

library(semPlot)

library(lavaan)

    • (2)練習1
      • 例1:線形回帰
        • 構造方程式(xがyを説明する場合)はy~x
        • 観測方程式(潜在要因fが観測量xを説明する場合)はy=~xとかく
        • 相関なしは潜在要因f1とf2が無相関はf1~0*f2とする

#通常の回帰

x<-rnorm(1000,mean=0,sd=2)

y<-3*x+rnorm(1000,mean=0,sd=0.1)

model1<-"y~x"

X<-cbind(x,y)

getCov(X)

#データ直接入力モード

fit<-sem(model=model1,data=X,estimator="ML")

#サマリー

summary(fit)

#サマリー、精度評価指標つき

summary(fit,fit.measure=TRUE)

#解析結果のグラフ化

semPaths(fit,whatLabels="stand",optimizeLatRes=T)

#解析結果のとりだし

parameterEstimates(object=fit)

#解析結果のとりだし標準化

standardizedSolution(object=fit)

      • 例3:二潜在要因を探す(体重(g)の列だけ除いて解析[相関係数1の要因があるとcov(X)が特異になり正定値行列でなく解析不明になるため])

#体重gの列を除く(体重(kg)と相関係数1の関係)

rr2b<-(rr2[,-which(colnames(rr2)=="体重(g)")])

#分散分散行列

r_d<-cov(rr2b)

r_d3<-r_d

#データのラベルを取得

n1<-labels(rr2b)2

#すべてのデータラベルの観測方程式の(y1=~長距離走; y1=~ 短距離走; y1~=ソフトボールなげなど)を作る。

model_test<-paste("y1=~",n1,sep="")

#すべてのデータラベルの観測方程式の(y2=~長距離走; y2=~ 短距離走; y2~=ソフトボールなげなど)を作る。

model_testb<-paste("y2=~",n1,sep="")

#潜在要因y1とy2は無相関という条件をいれる

model_testc<-'y1~0*y2'

#モデルをまとめる

model_test2<-paste(c(model_test,model_testb,model_testc),collapse="\n")

#データ数の取得

data_no<-length(rr2b[,1])

#SEMを実行(分散分散行列モード、データ数の指定が必要)

fit<-sem(model=model_test2,sample.cov=r_d3,estimator="ML",sample.nobs=data_no)

#semPaths(fit,whatLabels="stand",optimizeLatRes=T,nCharNodes=12)

#parameterEstimates(object=fit)

#標準化回帰係数の取り出し

standardizedSolution(object=fit)

    • 実際の解析
      • とりあえず、探索的因子分析(因子分析の結果)を利用して観測方程式を作ってみる。
      • 潜在要因をさらに朝ごはんや学校成績で説明するモデル(構造方程式

距離走<-(持久力),(瞬発力)

短距離走<-(瞬発力)

ソフトボール投げ<-(瞬発力)

朝ごはん->(持久力)

学校成績->(持久力)

運動習慣->(持久力)

朝ごはん->(瞬発力)

学校成績->(瞬発力)

運動習慣->(瞬発力)

朝ごはん->(筋力)

学校成績->(筋力)

運動習慣->(筋力)

*1:v3$loadings

2017-01-15

walfram alpha や mathmatica の整数条件の付け方の例

walfram alpha やmathmaticaの整数条件の付け方の例

スラッシュ「/;」とElement関数

Hypergeometric2F1[1-n+1/2,n+1/2,2-n+1/2,z] /; Element[n, Integers] && n > 0

dplyr ですべての列の平均を計算

dplyr ですべての列の平均を計算

dplyr::summarise_each(data,funs(mean),everything())

dplyr ですべての列の平均を計算

dplyr ですべての列の平均を計算

dplyr::summarise_each(data,funs(mean),everything())

2017-01-08

Rでキー列ごとに代表の値をえたデータフレームをえる方法。(Rでキーがユニークになるように集約、 グループ代表値以外データフレームから削除する方法)

Rでキー列ごとに代表の値をえたデータフレームをえる方法。ユニークになるように集約

つまり、グループ代表値以外データフレームから削除する方法は、distinctを用いる.

http://a-habakiri.hateblo.jp/entry/2016/11/29/215013

例:入力

a b c

1 1 3

1 1 4

1 2 5

2 1 6

2 1 7

2 1 8

3 1 9

4 2 10

aとbをキーとして代表値をして以下のテーブルがえたい

a b c

1 1 3

1 2 5

2 1 7

3 1 9

4 2 10

Rでは

library(dplyr)

a<-c(1,1,1,2,2,2,3,4)

b<-c(1,1,2,1,1,1,1,2,1)

c<-c(3,4,5,6,7,8,9,10)

dat<-data.frame(a,b,c)

dplyr::distinct(dat,a,b,.keep_all=TRUE)