RからLatexへ

メモ。今後、情報を更新していく。

私は完全にアンチ商用パッケージなので、できるだけビルゲイツにはお世話になりたくない。また、IBMにもお世話になりたくない*1。よって、分析はR、最近は文章をlatexで書くようになった。また、Rとlatexはとても互換性が良い。

分析結果の出力をR内のコードでlatex形式で出力できる。こうすれば、いちいち分析結果をコピペする必要がなく、ミスも減る。

よって、自分用に便利なコードの基本形をここにメモしていく（完全にpublication readyではないが、出来るだけ自動化する事を試みる）。分析の種類によっても載せたい統計量が異なるので今後更新していく。

＜記述統計量＞

#擬似データ作成
y<-rnorm(1000,mean=50,sd=10)
x1<-rnorm(1000,mean=50,sd=10)
x2<-ifelse(runif(1000)<=0.5,0,1)
data<-data.frame(cbind(y,x1,x2))

#ここから作成開始
num<-apply(data,2,length) #度数
mean<-round(apply(data,2,mean),3) #平均値
sd<-round(apply(data,2,sd),3) #標準誤差
min<-round(apply(data,2,min),3) #最小値
max<-round(apply(data,2,max),3) #最大値

smry_table<-cbind(num,min,mean,max,sd) #表形式に変換
colnames(smry_table)<-c("N","最小値","平均","最大値","標準偏差") #列名
rownames(smry_table)<-c("従属変数","独立変数","ダミー変数") #行名

library(Hmisc)
latex(smry_table, #当該テーブル
      file="sample.tex", #出力ファイル名
      title="", #1行1列目の文字
      caption="記述統計量" #題名（latexのcaption）
      )

＜回帰分析（パッケージ使用）＞*2

#擬似データ作成
y<-rnorm(1000,mean=50,sd=10)
x1<-rnorm(1000,mean=50,sd=10)
x2<-ifelse(runif(1000)<=0.5,0,1)
data<-data.frame(cbind(y,x1,x2))

#ここから作成開始
library(texreg)
texreg(lm(y~x1+x2,data=data)) #10%水準の有意は拾ってくれないので注意!

＜普通の回帰分析＞

   #擬似データ作成
y<-rnorm(1000,mean=50,sd=10)
x1<-rnorm(1000,mean=50,sd=10)
x2<-ifelse(runif(1000)<=0.5,0,1)
data<-data.frame(cbind(y,x1,x2))

#ここから作成開始
p_judge<-function(p){
  if(0.05<=p & p<0.1){rec<-"†"}
  if(0.01<=p & p<0.05){rec<-"*"}
  if(0.001<=p & p<0.01){rec<-"**"}
  if(p<0.001){rec<-"***"}
  if(0.1<=p){rec<-"n.s."}
  print(rec)
} #有意水準を判断する関数

rec<-summary(lm(y~x1+x2)) #分析結果収納
beta<-round(rec$coefficients[,1],3) #係数取り出し
sd.error<-round(rec$coefficients[,2],3) #標準誤差取り出し
p.value<-rec$coefficients[,4] #係数のp値取り出し
coef.sig<-numeric(length(p.value))
for(i in 1:length(p.value)){
  coef.sig[i]<-p_judge(p.value[i])
} #係数の有意水準
adj.r.squared<-round(rec$adj.r.squared,3) #自由度調整済みr2
f.test<-df(x=rec$fstatistic[1], df1=rec$fstatistic[2], df2=rec$fstatistic[3]) #f検定のp値
f.test.sig<-p_judge(f.test) #f検定の有意水準

table<-cbind(beta,sd.error,coef.sig) #表形式に変換
colnames(table)<-c("回帰係数","有意水準","標準誤差") #列
table<-rbind(table, adj.r.squared,f.test.sig)

library(Hmisc)
latex(table, #当該テーブル
      file="sample.tex", #出力ファイル名
      title="", #1行1列目の文字
      caption="回帰分析の結果" #題名（latexのcaption）
)

＜glm系全て＞ 

y<-ifelse(runif(1000,min=0,max=1)<0.5,0,1)
x1<-rnorm(1000,mean=50,sd=10)
x2<-rnorm(1000,mean=50,sd=10)

#ここから作成開始
p_judge<-function(p){
  if(0.05<=p & p<0.1){rec<-"†"}
  if(0.01<=p & p<0.05){rec<-"*"}
  if(0.001<=p & p<0.01){rec<-"**"}
  if(p<0.001){rec<-"***"}
  if(0.1<=p){rec<-"n.s."}
  print(rec)
} #有意水準を判断する関数

fit<-glm(y~x1+x2,family="binomial")

rec<-summary(fit) #分析結果収納
beta<-round(rec$coefficients[,1],3) #係数取り出し
sd.error<-round(rec$coefficients[,2],3) #標準誤差取り出し
p.value<-rec$coefficients[,4] #係数のp値取り出し
coef.sig<-numeric(length(p.value))
for(i in 1:length(p.value)){
  coef.sig[i]<-p_judge(p.value[i])
} #係数の有意水準
AIC<-round(rec$aic,3) #AIC取り出し
BIC<-round(BIC(fit),3) #BIC取り出し
library(BaylorEdPsych) 
pseudo.R2<-round(PseudoR2(fit)[1],3) #擬似決定係数取り出し

table<-cbind(beta,sd.error,coef.sig) #表形式に変換
colnames(table)<-c("回帰係数","標準誤差","有意水準") 
table<-rbind(table, AIC, BIC, pseudo.R2)

library(Hmisc)
latex(table, #当該テーブル
      file="sample.tex", #出力ファイル名
      title="", #1行1列目の文字
      caption="回帰分析の結果" #題名（latexのcaption）
)

＜t検定（等分散仮定）＞

####t検定####
#擬似データ作成
y<-rnorm(1000,mean=50,sd=10)
x<-ifelse(runif(1000)<=0.5,0,1)
data<-data.frame(cbind(y,x))

#ここから作成開始
group1<-subset(data,data$x==1)
group0<-subset(data,data$x==0)

result<-t.test(group1$y,group0$y, var.equal = T)

p_value<-round(result$p.value,3)
t_value<-round(result$statistic,3)
group1_mean<-round(mean(group1$y),3)
group0_mean<-round(mean(group0$y),3)

p_judge<-function(p){
  if(0.05<=p & p<0.1){rec<-"†"}
  if(0.01<=p & p<0.05){rec<-"*"}
  if(0.001<=p & p<0.01){rec<-"**"}
  if(p<0.001){rec<-"***"}
  if(0.1<=p){rec<-"n.s."}
  print(rec)
}

sig<-p_judge(p_value)

t_table<-cbind(group1_mean, group0_mean,t_value,p_value,sig)
colnames(t_table)<-c("グループ１の平均値","グループ0の平均値","t値","p値","有意水準") #列名
rownames(t_table)<-c(" ") #行名

library(Hmisc)
latex(t_table, #当該テーブル
      file="sample.tex", #出力ファイル名
      title="", #1行1列目の文字
      caption="t検定の結果" #題名（latexのcaption）
)

*1:最近、研究室で予算がないとか言っているが、まずは商用パッケージを最小限にしたらどうかと思う。これいくら言っても反感を買うので、もう言わないようにしている。100歩譲ってwordは理解できるが、spssはもうやめたらどうかと思う。プルダウンで分析したいのなら最近はHADとかあるし。

*2:stargazerなどもあるらしいが、とりあえずtexregが簡素で応用しやすそう。また一般化線形モデルのglm()や順序ロジットのclm()（ordinalパッケージ）にも対応している。またマルチレベル分析のlmer()（lme4パッケージ）やマルチレベルの一般化線形モデルのglmer()（lme4パッケージ）、マルチレベルの順序ロジスティック回帰分析のclmm()（ordinalパッケージ）にも対応している。