杉本隆さん（おすぎ）から提供いただいた、前回のRを動かした回のメモをアップします！
おすぎ、ありがとう！！
皆様、書式が少しずれてしまっているのですが、まずはアップしてしまいます。
更新：2011/06/09
【コンセプト】
□データ解析の実務者の視点からRを使用
□データを取り扱う際に遭遇する困難を明らかに
□実例プログラムは必ず入れる
□安定第一という観点からR本体で全て進める（コマンダーは参照程度）
0. プログラムを書き始めるために
Rを起動して最初に出てくるR Consoleは，結果の出力と一緒になっていて，プログラムを加工しにくい。
メニューの[ファイル]→[新しいスクリプト]でRエディタを起動する。
エディタに書いたプログラムを実行するときは，選択して[Ctrl+r]をするか，R Consoleにコピペして実行
1. データの読み込み
データの種類には色々あるが、現実場面でよく使われそうな形式として、Excelに入力したデータ（csvファイルに変換したもの）やSPSSのデータファイル（savファイル）を取り上げる。まず，csvファイルの読み込みを解説する。csvファイルを使うのは，RでExcelファイル（xlsファイル）をそのまま読もうとすると大変な理由による。
まず，読み込むときにフォルダのパス名を教えてあげなければいけないのだが，普段はこれは面倒すぎる
ので，省略するために，「作業ディレクトリ」という，データの置いてあるフォルダの指定を行う
作業ディレクトリを変更する場合は
メニューの[ファイル]→[ディレクトリの変更]でフォルダを指定。デスクトップにフォルダを作ると分かりやすいかも。
Rは，データを「オブジェクト」として認識する。基本的には，<-の矢印が向けられた先にあるものが
オブジェクトとなる。名前は自由。
そして，今度はcsvファイルをRに読み込む。これからデータ加工を始めるので，順番に作っていくことを考えて，オブジェクト名を「a」とした。
作業ディレクトリを指定しているので，"jgss_mini.csv"と，データファイル名だけでOK
a <- read.csv("jgss_mini.csv", header=T) # header=Tは一行目が変数名という意味
ファイル名打つのが面倒ならば
a <- read.csv(file.choose(), header=T)
読み込んだデータセット確認のため
view(a)→デフォルトでは×
head(a) # 最初の6人表示
fix(a) or edit(a)
オブジェクト確認
ls()
オブジェクト削除するときは
rm(a)
全部削除は
rm(list=ls())
どうせ1つのデータセットで集中的に作業するんだから，attachでいいんじゃないか。
-----------------------------------------------
とりあえず，度数分布でも確認してみたい
table(a$DOPET)
%を出すには
prop.table(table(a$DOPET))
クロス
table(a$DOPET,a$SEXA)
library("Hmisc")
describe(a)
describe(a$DOPET)→n，欠損値，平均，度数分布まで出る。factorだろうがnon-factorだろうが関係ないよ
library("gmodels")
Frequencies and Crosstabs
http://www.statmethods.net/stats/frequencies.html
--------------------------------------------------
2. データマネジメント
(1)数量型変数を因子型変数に変換，またその逆
全て数量変数として読み込まれているので，カテゴリとして扱いたい変数を因子型に変換する。Rはこれらの区別を厳密につけたがるため
a$CASEID <- as.factor(a$CASEID)
a$SEXA <- as.factor(a$SEXA)
a$MARC <- as.factor(a$MARC)
a$XJOB1WK <- as.factor(a$XJOB1WK)
a$DOPET <- as.factor(a$DOPET)
a$DODOGE <- as.factor(a$DODOGE)
a$DODOGI <- as.factor(a$DODOGI)
a$DOCAT <- as.factor(a$DOCAT)
a$DOMAMMAL <- as.factor(a$DOMAMMAL)
a$DOBIRDA <- as.factor(a$DOBIRDA)
a$DOFISHA <- as.factor(a$DOFISHA)
a$DOREPTL <- as.factor(a$DOREPTL)
a$DOINSECT <- as.factor(a$DOINSECT)
a$DOOTHER <- as.factor(a$DOOTHER)
a$DOPNOMK <- as.factor(a$DOPNOMK)
summary(a)
で確認
以上は1つずつ変換したが，↓の3行を使えば，それだけでできる。
作業内容は，オブジェクト該当変数をfactorに変換した「af」を作り，aとmergeして，新データセット「b」を作成
af <- sapply(data.frame(a[,1:2], a[,4:6], a[,8:19]), factor)
af <- data.frame(af)
b <- merge(af, a, by=names(af), all=T) #aとafの順番を間違えないように
summary(b)
で確認
因子から数値，NAが入っている数値データはそれをやる必要性あり
(2)欠損値の定義
(3)欠損値の除去
(4)連続変数をカテゴリ変数に（例：年齢層作成）
(5)逆転項目の定義
--------------------------------
[応用]
(1)
SPSSのデータファイルを読み込んで、汎用性の高いcsvに変換して保存し，さらにそれを読み込むという作業
まず、他のデータ形式を読み込むための準備として、foreignパッケージを読み込む
library(foreign)
そして、データファイルの置いてあるフォルダを指定してデータセット（jgss_mini.sav）を読み込む
ここでは「jgss」という名前のオブジェクトで
jgss <- read.spss("データファイルのフォルダまでのパス/jgss_mini.sav", to.data.frame=T, use.value.labels=F)
しかし，パス指定は最初はハードルが高いので，以下の作業ディレクトリの指定の方法を薦める
*[作業ディレクトリの場所指定が面倒くさい場合]
作業ディレクトリの場所にデータファイルがあれば，フォルダのパス名を省略できる。
現在の場所は
getwd()
で確認できる
プログラムでは，
setwd("ここにフォルダのパス/ファイル名")
であり，フォルダのパスは，フォルダのアドレス欄からコピペできるが，
Windowsでは，フォルダの区切り記号が「\」なため，これを「/」に変換する必要
があり，少し面倒くさい。
これで、場所指定が省略されてファイル名だけで作業できるようになる。
jgss <- read.spss("jgss_mini.sav", to.data.frame=T, use.value.labels=F)
--------------
次に，読み込んだsavファイルをcsvファイルに変換して保存する
write.table(jgss, "jgss_mini.csv", sep=",", row.names=FALSE)
あとはcsvの読み込みと同じ
20110610 mixed model multilevel model
> results1 <- lmer(mathach ~ 1 + (1| schid), data = a )
> #これはnull modelを表現している！
> ＃スクールによって切片が違うかを見ている
>
> summary(results1)
Linear mixed model fit by REML
Formula: mathach ~ 1 + (1 | schid)
Data: a
AIC BIC logLik deviance REMLdev
47123 47143 -23558 47116 47117
Random effects:
Groups Name Variance Std.Dev.
schid (Intercept) 8.614 2.9350
Residual 39.148 6.2569
Number of obs: 7185, groups: schid, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.6370 0.2443 51.72
> #deviance値でモデルの適合度を判断するのが通例らしい
> ＃Interclass correlation coefficient (ICC)＝8.614/8.614+39.148である。σ＾2＝39.148。
> ＃ICC＝0.18はだいぶ高いらしい。学校によって、18%の効果がある→mixed modelを使用する価値がある
> ＃学校の数はある程度（20とか？）ないと、信頼性が怪しくなる（３校とかでは怪しいのではないか・・）
>
> results2 <- lmer(mathach ~ 1 + meanses + schtype + (1 | schid), data = a)
> summary(results2)
Linear mixed model fit by REML
Formula: mathach ~ 1 + meanses + schtype + (1 | schid)
Data: a
AIC BIC logLik deviance REMLdev
46956 46991 -23473 46944 46946
Random effects:
Groups Name Variance Std.Dev.
schid (Intercept) 2.3139 1.5211
Residual 39.1614 6.2579
Number of obs: 7185, groups: schid, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.0963 0.1986 60.90
meanses 5.3329 0.3686 14.47
schtype 1.2253 0.3058 4.01
Correlation of Fixed Effects:
(Intr) meanss
meanses 0.245
schtype -0.698 -0.356
> #devianceが減っているので、より適切なモデルであることがわかる
> ＃ランダム切片も、8.614から2.31に減っているのでより説明できている
> ＃meansesが１ポイント上がると、5点くらい数学の点数が上がる
> ＃lme4の作者が、自由度が複雑なので、Ｐ値は出さないように設計した
> results3 <- lmer(mathach ~ meanses + schtype + ses + meanses*ses + schtype*ses + (1 + ses | schid), data = a )
> summary(results3)
Linear mixed model fit by REML
Formula: mathach ~ meanses + schtype + ses + meanses * ses + schtype * ses + (1 + ses | schid)
Data: a
AIC BIC logLik deviance REMLdev
46525 46594 -23253 46498 46505
Random effects:
Groups Name Variance Std.Dev. Corr
schid (Intercept) 2.407424 1.55159
ses 0.014619 0.12091 1.000
Residual 36.758381 6.06287
Number of obs: 7185, groups: schid, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.0948 0.2028 59.65
meanses 3.3202 0.3885 8.55
schtype 1.1949 0.3079 3.88
ses 2.9052 0.1483 19.59
meanses:ses 0.8462 0.2718 3.11
schtype:ses -1.5781 0.2245 -7.03
Correlation of Fixed Effects:
(Intr) meanss schtyp ses mnss:s
meanses 0.213
schtype -0.676 -0.344
ses 0.076 -0.146 -0.064
meanses:ses -0.143 0.178 -0.081 0.278
schtype:ses -0.062 -0.080 0.093 -0.679 -0.356
＃Rstudioをインストールすると便利！
