gikoha’s blog

個人的メモがわり

統計解析

データを増やしたので解析に入るが、R + mysqlを使った方が解析は速い感じがする。
特にNA、欠損データの扱いがやっぱりうまい。
ただ、データベースでint/float型にしてしまうと、空欄やNAと入力したところが0になってしまう。
データベースを作るときにはVARCHARにしておいて、数値、空欄、"NA"含めてみんな文字としてデータベースに代入。

library(RMySQL)

drv<-dbDriver("MySQL")
conn<-dbConnect(drv,user="XXX", password="XXX", dbname="XXXdb")
d<-dbReadTable(conn,"XXXtable")
dbDisconnect(conn)

こうやって読み込んだ後、

# NA、空欄を NAとし、それ以外を数値へ変換する

d$BH = as.numeric(d$BH)
d$BW = as.numeric(d$BW)

数値が必要なrowだけ数値に変換

でもグラフはgplot2とかでチマチマ書いていくよりprismに書き戻したほうが絶対きれいですな

だからsummarySE

summaryPrism <- function(data=NULL) {
	size= dim(data)[1]
	vec <- c()
	for(i in 1:size){
		vec <- append(vec, data[i,3])		# value
		vec <- append(vec, data[i,4])		# sd
		vec <- append(vec, data[i,2])		# N
	}
	return (vec)
}

こいつを使って

> dfwc<-summarySE(d,measurevar="Lhratio", na.rm=TRUE, groupvars="ACS");
 summaryPrism(dfwc);
 t.test(Lhratio~ACS, data=d)
[1]  2.0498168  0.8261261 80.0000000  2.6501885  0.9033336 46.0000000

	Welch Two Sample t-test

data:  Lhratio by ACS 
t = -3.7041, df = 87.201, p-value = 0.0003715
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.9225159 -0.2782275 
sample estimates:
mean in group 0 mean in group 1 
       2.049817        2.650188 

裏でPrismからGrouped Tableを作成する

Rに戻って上の [1] 以下の「2.0498168 0.8261261 80.0000000 2.6501885 0.9033336 46.0000000」を
選択コピーして Prismの中のテーブルにペースト

グラフタイプはColumn mean ,error barに変更するともうグラフが書けている
pとかタイトルとかを適当に手直しすればあっというまにスライド/論文品質のグラフが書けました

ちなみにこれと似たことをgplot2でやろうとすると..

dfwc<-summarySE(d,measurevar="Lhratio", na.rm=TRUE, groupvars="ACS")
p<-ggplot(dfwc,aes(y=Lhratio,x=ACS))
p+geom_point()
 +geom_errorbar(aes(max=Lhratio+se,min=Lhratio-se,width=0.2))
    +scale_x_continuous(breaks=c(0,1),labels=c("ACS(-)","ACS(+)"))+ylab("LDL-C/HDL-C ratio")

これだけのコードを書いて、

このクオリティではダメだね

表もちまちまlatexでやっている人がいるけど、そんなんやりたくないわー
なんとかならんかな