R
# StataのデータとSPSSのデータ変換
# Stataのデータ（dta）をSPSSのデータ（sav）にRのhavenパッケージで変換
library(haven)
d <- read_dta("file_name.dta")
write_sav(d, "file_name.sav")
# SPSSのデータ（sav）をStataのデータ（dta）にRのhavenパッケージで変換
d <- read_sav("file_name.sav")
write_dt(d, "file_name.dta")
#クロス表の作成がめんどくさいのでリスト
#暇があれば多重クロス表まで
#もっと暇があればログリニアまで
library(magrittr)
cro <- function(Row,Column){
if(is.table(Row)) {
d <- data.frame(Row)
d <- d[rep(row.names(d), d$Freq), 1:2]
Row <- d[,1]
Column <- d[,2]
}
tbl <- table(Row,Column)
tblsum <- tbl %>% addmargins()
rowtbl <- tbl %>% addmargins(1) %>% prop.table(1) %>% addmargins(2) %>% round(3)*100
coltbl <- tbl %>% addmargins(2) %>% prop.table(2) %>% addmargins(1) %>% round(3)*100
chisq <- chisq.test(tbl)
V <- sqrt(chisq$statistic/sum(tbl)/(min(dim(tbl))-1))[[1]] %>% round(3)
list(N = sum(tbl),Table = tblsum, Row_Percent = rowtbl, Column_Percent = coltbl, Chi_squared_Test = chisq, V = V)
}
cro(occupationalStatus)
#期待度数
tbl <- matrix(c(1,2,6,3,6,7,5,2,1,4,5,6),nrow=3)
tbl
#Efre <- outer(rowSums(tbl), colSums(tbl), "*")/sum(tbl)
#Efre
Efre <- matrix(NA,nrow=nrow(tbl),ncol=ncol(tbl))
Efre
for (j in 1:ncol(tbl)){
for (i in 1:nrow(tbl)){
Efre[i,j] <- (sum(tbl[i,])*sum(tbl[,j]))/sum(tbl)
}
}
dimnames(Efre) <- dimnames(tbl)
Efre
Chisq <- sum((tbl-Efre)^2/Efre)
df <- (nrow(tbl)-1)*(ncol(tbl)-1)
p <- pchisq(Chisq, df, lower.tail = FALSE)
names(Chisq) <- "Chi-square statistic"
names(df) <- "Degree of freedom"
names(p) <- "p-value"
list(Obs=tbl,Exp=round(Efre,1),Chisq=Chisq,df=df,p=p)
chisq.test(tbl)
2014年9月4日（木）・5日（金）10:00～17:15から，「2014年 計量分析セミナー」「2014年度 第9回 ICPSR国内利用協議会・統計セミナー」
で，「Rによる二次分析入門」の講師を担当いたします．
http://ssjda.iss.u-tokyo.ac.jp/351010/pdf/R2014.pdf
SSJDAに寄託されている社会調査データ（SPSS形式）を，Rで読み込み，変数を加工し，分析を行うという
一通りのプロセスを，複数の多変量解析の紹介を交えながら実践します．
Rに関するメモ
【SPSSデータの保存】
SAVE TRANSLATE OUTFILE='D:\R\file.dat'
/TYPE=TAB /MAP /REPLACE /FIELDNAME
/CELLS=VALUES
KEEP= varA .
【Rでデータを読み込む】
x <- read.table("file.dat", header=T, sep="\t")
【ケース別に分析】
result <- by(x, x$q01_1, with,
lm(rHgaku ~ q01_2a + EGP6ana + edu_4 + loghinc)
)
summary(result[[1]])
summary(result[[2]])
【jitterをつかって，散布図を少し見やすくする】
#x座標もy座標もずらす
plot(jitter(x),jitter(y))
#x座標のみずらす
plot(jitter(x),y)
#y座標のみずらす
plot(x,jitter(y))
#たくさんずらす
plot(jitter(x,5),jitter(y,5))
【記述的分析】
#たとえば次のようなデータyがあるとする．従属変数はLifeSである．その他の変数は独立変数である．
#従属変数は一番前に，独立変数は因果の逆の順に並べるとよい．
y <- data.frame(LifeS,Hincome,OccPres,Eduy,Age)
#散布図，回帰直線，相関係数を書きたい場合，pairs関数を利用する．
#panel.lsfit2のところでjitterを使うと図がみやすく（単調で）なくなる．
panel.lsfit2 <- function(x,y,...) {
f <- lsfit(x,y)$coef
xx<- c ( min (x) , max(x) )
yy <- f["X"] * xx+ f["Intercept"]
lines(xx,yy, lwd=3)
points(jitter(x),jitter(y))
}
panel.cor2 <- function(x, y, digits=2, prefix="", cex.cor)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- round((cor(x, y)),2)
txt <- format(c(r, 0.123456789), digits=digits)[1]
txt <- paste(prefix, txt, sep="")
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = 2)
}
pairs(y, upper.panel=panel.lsfit2, lower.panel=panel.cor2)
#出来上がったのが以下の図である．
【3変数の関連を同時に表示】
太郎丸先生のブログのエントリをもとに作成．
http://sociology.jugem.jp/?day=20091208
用いるのはggplot2．http://had.co.nz/ggplot2/
#ライブラリの読み込みとデータセットの作成
library(ggplot2)
X <- (1:20)/2
Y <- X+rnorm(20, sd=2)
Z <- X+0.5*Y+rnorm(20)
Z <- round((Z - min(Z)+1)/4,1)
D<-data.frame(X,Y,Z)
summary(D)
#色で第3変数の大きさを表現
qplot(X, Y, data=D, colour=Z) +
scale_colour_gradient(low="pink", high="red",limits=c(min(Z),max(Z)))
#白黒で第3変数の大きさを表現
qplot(X, Y, data=D, colour=Z) +
scale_colour_gradient(low="grey80", high="black",limits=c(min(Z),max(Z)))
#数値と白黒で第3変数の大きさを表現
qplot(X, Y, data=D, colour=Z,label=factor(Z), geom="blank") +
geom_text()+ scale_colour_gradient(low="grey75", high="black",limits=c(min(Z),max(Z)))
最後のものが，情報量が多くていいか．
ちなみにggplot2で作成した図の背景は以下のようにすることでいじることができる．
old_theme<- theme_update(
panel.background = theme_rect(fill = "grey85"),
axis.text.x = theme_text(colour = "black",size=12, vjust = 1),
axis.text.y = theme_text(colour = "black",size=12, hjust = 1),
axis.title.x = theme_text(colour = "black", face = "bold",size=14, vjust = 1, hjust = 0.4),
axis.title.y = theme_text(colour = "black", face = "bold",size=14, angle = 90,vjust = -0.2 ))
出来上がり．
