心理系の大学生が、Rをそこそこ使えるようになるサイトを目指して・・・

#######################1要因分散分析############################## anova(aov(従属変数~独立変数)) anova(aov(従属変数~独立変数,data=データ名)) oneway.test(従属変数~独立変数,data=データ名,var.equal=T) #var.equal=T　は各データの等分分散性の仮定を行うための設定 #######################2要因分散分析(繰り返しのない）############################## summary(aov(従属変数~独立変数1+独立変数2,data=データ名))

#######################重回帰分析############################## data_fin<-c(factor1,factor2,factor3,factor4)   #使いたい変数を1つのデータにまとめる class(data_fin)   #上でまとめたデータの型を確かめる data_fin<-data.frame(data_fin)   #データの型を”data.frame”型に変更する data_fin.sd<-scale(data_fin)   #データ内の変数を標準化する class(data2_fin.sd)   #また、データの型を確認する data_fin.sd<-data.frame(data_fin.sd)   #"data.frame"型に変更する summary(lm(formula=factor1~factor2+factor3+factor4,data=data_fin.sd)) #summary(lm(従属変数~独立変数1+独立変数2+独立変数3,data=データセット名)) #結果の見方 #帰無仮説は「係数＝0」 #AICは値が小さい方がモデルの当てはまりの良さを示す #Multiple R-squaredが決定係数 # library(DAAG) result<-lm(formula=factor1~factor2+factor3+factor4,data=data_fin.sd) vif(result) #多重共線性の確認 #vif関数を使うには、”DAAG”パッケージが必要です。  #VIF>10だと多重共線性に問題あり。<2なら問題なし。 ##############ステップワイズ変数選択による重回帰分析###################### sreg("f2s","f3s","f4s","f5s","a1", stepwise=TRUE, P.in=0.05, P.out=0.05, pre

#######anova kun########## source("anovakun_480.txt") #anovakunを読み込む #anovakunのデータを指定したディレクトリに入れておく必要がある。 #anovakunのファイル名は、versionにより異なるようなので、自分の持っているanovakunが入ったファイル名を確認してください。 aov_1<-data[,c("要因A","要因B","従属変数")] #分析専用のデータセットを作成、データセットの最後が従属変数になる #注）分析ごとにデータセットを作らなければならない anovakun(データセット,"要因計画",水準数,holm=T,peta=T) #anovakun(データフレーム名,要因により変わる(*1),水準数,多重比較(ホルム法)の有無,効果量の指標の設定) #(*1)：2要因だとAB,3要因だとABC #被験者間←s→被験者内 （例 １要因被験者内：sA ２要因被験者間：ABs ３要因混合計画：ABsC #要因計画の書き方-> 被験者間要因は「s」の左側，被験者内要因は「s」の右側に配置します。 #例えば，２要因の被験者間計画は「ABs」，３要因の被験者内計画は「sABC」 #水準数の書き方はデータセットの要因の順番に対応する #要因の指定の関係で、被験者内要因を最後にデータセットを作成する必要がある。 ##以下、オプションについて## #holm・・・「holm = T」とすることで，Holmの方法（Sequentially Rejective Bonferroni Procedure）を指定できます。Shafferの方法よりも若干検出力が落ちますが，使用に際しての制約がより少なく，広い場面で適用できる方法です。 #hc・・・「hc = T」とすることで，Holland-Copenhaverの方法（Improved Sequentially Rejective Bonferroni Test Procedure）を指定できます。Shafferの方法よりも若

########カテゴライズ###################### #3郡に分ける。平均値±1/2SDで区切る data_n　<-　data2_uwaki_n0 describe(data_n) result<-describe(data$ data$a01 ) low<-result$mean-0.5*result$sd high<-result$mean+0.5*result$sd data$a01_cate[data$a01<=low]<-1 data$a01_cate[ data$a01 >low& data$a01 <high]<-2 data$a01_cate[ data$a01 >=high]<-3 ##############cut関数（カテゴリライズ） cut(データ名,breaks=c(0,80,100), labels=c("不合格","合格"), right = FALSE, include.lowest = TRUE) #breaks=c(0,80,100)←区切り位置の指定（ラベル数よりも1つ多くなる） #right ←デフォルトは(TRUE) #0は入らず80まで（以下）、80は入らず100まで #FALSEだと、0から80未満、80以上100未満 #include.lowest = TRUE　←　right = FALSEの場合は100がカテゴリに含まれないため、100（上限）をカテゴリに含む設定 私は個人的にcut関数が好きですが、”以上”や”未満”の設定が面倒なのと、 柔軟性にかける部分もあると思うので（私の勉強不足かもしれませんが...） 上段に記載した方法が、アンパイです。

#### t 検定 t.test(a1~data2$sex) #t.test(従属変数~独立変数) #### t 検定Part2（繰り返しpart1） t.test2<-function(x,y1,y2,z){   for(a in y1:y2){     cat("従属変数")     print(names(x[a]))     cat("独立変数")     print(names(x[z]))     print(t.test(x[,a]~x[,z]))   } } t.test2(data2,26,62,2) #### t 検定Part3（繰り返しpart2（未検証）） t.test3<-function(x,z,y1,y2){   for(a in y1:y2){     cat("従属変数")     print(names(x[z]))     cat("独立変数")     print(names(x[a]))     print(t.test(x[,z]~x[,a]))   } } t.test3(data,,,) ###ちょっとしたコツ t 検定をやる時に、ついついデータセットを別にしたくなる人がいるようですが、そんな必要はありません。 独立変数となる項目がしっかりと元のデータセットにあれば、問題ありません。 むしろ、データセットは乱立すると、よくないと思います。欠損処理や結合時に行がズレるなどの問題の原因になります。 独立変数はt検定ですので、２値である必要があります。２値であればfactor型出なくてもおそらく分析できます。 ###エラー 独立変数が２値でない場合に、エラーが出ます。３値以上のデータである、もしくはデータが２値未満であることが考えられます。table関数を使って、どのようなデータになっているか確認するといいでしょう。

###################因子分析編######################### VSS.scree(xx) #スクリープロット作成 x.cor<-cor(xx,use="complete.obs") #相関マトリックスを準備して #[use="complete.obs"]=>欠損を含むサンプルをすべて除く #[use="pairwise.complete.obs"]=>2変数の組み合わせごとに欠損値があるサンプルを除く（ペアワイズ） eigen(x.cor) #因子ごとの固有値$values fa.parallel(x.cor) #平行分析 d1<-fa(x.cor,nfactors=5,fm="ml",rotate="promax") #nfactors=抽出する因子数 #fm=因子抽出方法(pa=主因子法,ml=最尤法,gla=一般化最小二乗法) #rotate=回転の種類(promax=プロマックス,varimax=バリマックス,none=なし) print(d1,sort=TRUE,digit=3) #sort=TRUE並び替えの有無 #digit=負荷量（パターン）の表示桁数 #h2共通性 #u2独自性 #SS loadings負荷量の二乗和 #Proportion Var寄与率 #Cumulative Var累積寄与率 ################因子分析編Part2（結果を負荷量の低いものを非表示にする）##################### #resultというデータセットをつくり、その中に因子分析結果をぶちこむ result <- fa(xa10,nfactors=5,fm="ml",rotate="promax") #因子分析(データセット名,因子数=4,抽出法＝ml[mlは最尤法。],回転=プロマックス) fa(xa10,nfactors=5,fm="ml",rotate="p

##############α係数################### F1 <- data[,c(“a1”,”a2”,”a3”,”a4”)] alpha(因子) #raw_alpha　α係数 #std.alpha　標準化されたα係数 #G6(smc)　ガットマンのラムダ6 #average_r　項目間相関の平均値 #mean　平均 #sd　標準偏差 #############因子得点################ attach(data) #attach(データ名)　←　データの中の項目を検索可能にする("data$"を項目の前につける必要がなくなる) factor_score <- (項目A+項目B+項目C)/3 #全ての因子項目を足し、項目数で割る 項目数を割る処理は、してもしなくてもどちらのパターンもあります。 相関などを見る際は、項目数で除算しなくても行えます。 平均値の差を見る検定（t検定や分散分析）では、項目数で除算しないと駄目です。

cor(項目名,項目名) #.00~.20ほとんど相関がない #.20~.40弱い相関がある #.40~.70比較的強い相関がある #.70~1.0強い相関がある cor.test(項目名,項目名) #無相関検定 data_a<-data2[,3:25]   #data2の3列目から25列目をdata_aに代入 soukan<-cor(data_a)   #cor(data_a)←相関行列を出す   #それを、soukanに代入 write.csv(soukan,"tmpcor1.csv")   #soukanを"tmpcor1.csv"の名前で出力   #第一引数はR内でのデータ名、第二引数は書き出すファイル名 round(cor(データ名,use="pairwise.complete.obs"),3)   #round=>小数点以下の桁数を指定   #pairwise.obs=>ペアで使うので、最大限にデータを使える   #complete.obs=>欠損のある行は使わない

install.packages("scatterplot3d") library("scatterplot3d") trees <- data[,c("age_kc","age","pre_f1s")] scatterplot3d(trees) #3次元の動かない図が出力できる install.packages("rgl") library("rgl") plot3d(trees$age_kc, trees$age, trees$pre_f1s) #3次元で出力後、グリグリ動かせる writeWebGL(width=50, height=50)

###############知ってると便利################### #Rを終了 q() #作業スペースにあるオブジェクト名をみる ls() #a&b&cのオブジェクトを削除する rm(a,b,c) remove(a,b,c) par ( mfrow=c ( 2 , 2 ))      #2*2で図を配置してくれる。 hist(変数名)     #ヒストグラム pie(table(因子形式の変数名))     #円グラフ boxplot(変数名)     #箱ひげ図 boxplot(変数名~条件変数名)     #層別の箱ひげ図 plot(変数名1,変数名2)     #散布図 pairs(複数の変数名の指定)     #散布図行列 #クロス集計表 table( データフレーム名 $ 変数名m, データフレーム名 $ 変数名n ) table( データフレーム名 [ c ( 列番号m, 列番号n ) ] ) #x内にある変数名を直接取り扱えるように、検索パスに追加する attach(a) #xを検索パスから削除する detach(x) #xにある変数aの型を確認する class(x$a) #データセットxに、作られた変数new1とnew2を横につなげてひとつにする x<-cbind(x,new1,new2) #条件式にあうケースのみを取り出す xnew<-subset(x,条件式) #（すべてのケースの）変数名1、変数名2、変数名3のデータを指定する x[,c("変数名1","変数名2","変数名3")] #変数名を0のみ選ぶ 変数名==0 #変数名を男性、のみを選ぶ 変数名=="男性" #"ctrl+F"で検索ができる #"ctrl+Alt+b"で全部実行

##############基本作業編################## setwd("ディレクトリ")   #ディレクトリの設定 data <- read.csv("ファイル名.csv")   #ファイルの読み込み   #注意）スラッシュの向きに注意 data2 <- na.omit(data)   #dataから欠損値のあるデータを行ごと削除して、data2 に代入   #データによっては、欠損処理によってサンプル数が減り過ぎるので、   #欠損処理のタイミングと方法には注意 count<-function(x) {sum(!is.na(x))}   #データの欠損以外の個数を数える式 count(data$rikai) by(data2$sex,data2$sex,count)   #by(データ名,群名,何をするか) describeBy(データ名,group = data$sex)   #"group = グループ分けしたい変数名"   #グループ分けして、基礎統計量を出す describe(data2)   #data2の基礎統計量を出す ncol(x) #行数 nrow(x) #列数 data2_m<-subset(data2,sex==1) data2_f<-subset(data2,sex==2) #data2からsex1.2ごとにそれぞれdata2_mとdata2_fに代入 #"&"でAND、"|"でOR table(data$gra) #度数分布 hist(data$gra) #ヒストグラム

############## 共分散構造分析( lavaan ) ############################ コード概要 モデル名 <- # 測定変数の指定( =~ を使う) 潜在変数名 1=~ 観測変数 1+ 観測変数 2+ ・・・・・ 潜在変数名 2=~ 観測変数 1+ 観測変数 2+ ・・・・・ # 回帰式( ~ を使う) 従属変数名 1~ 説明変数 1+ 説明変数 2+ ・・・ 従属変数名 2~ 説明変数 1+ 説明変数 2+ ・・・ # 分散、共分散( ~~ を使う) 変数名 i~~ 変数名 i # 分散 変数名 j~~ 変数名 k # 共分散 モデル部分はシングルカンマ( ' )でくくる すべての変数の分散を推定する場合は分散の指定は省力できる パラメタ値の推定 lavaan オブジェクト名 <-sem (モデル名, data= データ名） または lavaan オブジェクト名 <-lavaan (モデル名, data= データ名, model.type="sem" ,    auto.var=TRUE , auto.fix.first=TRUE , auto.cov.lv.x=FALSE , auto.cov.y=FALSE ) #auto.var=TRUE      : すべての変数の分散または残差分散を推定する #auto.fix.first=TRUE: 非標準化解において測定変数のうち最初の 1 つのパス係数を 1 に固定する #auto.cov.lv.x=FALSE: 潜在変数間の共分散を自動的には推定しない #auto.cov.y=FALSE   : 従属変数の残差間の共分散を自動的には推定しない #int.ov.free ： 観測変数の切片を自由にするか（ FALSE なら 0 に固定） #int.lv.free ： 潜在変数の切片を自由にするか( FALSE なら 0 に固定） summary ( lavaan オブジェクト名, fit.measures=TRUE , standardized=TRUE ) # 出力オプション fit.measur