2016/07/22

2020/05/04

実行するだけでR言語入門出来るプログラム

R言語入門

ライター:

おまけとして、入門コーナーの内容を一気に勉強出来るプログラムを作りました。プログラムを1行ずつ実行するだけで、R言語の基本的な学習が済むというような内容になっています。時間がない人や手っ取り早く勉強を入門を済ませたい人におすすめです。

作業ディレクトリに以下の2つのファイルがあれば、プログラムを流すことが出来ます。

sample-data

cox_sample_data

解析用に作ったサンプルファイルですので、ない人はダウンロードしておきましょう。

※R言語入門のトップページはこちら

入門用プログラム

プログラム

###ベクトル・行列、データの読み込み、データ整理基礎###
#ベクトルの作成
x <- c(1,2,3,4,5)
x
y <- c(1:5,3:1)
y
z <- c(rep(3,4),rep(c(1,5,10),c(2,3,4)))
z
a <- c(“A”,”B”,”C”)
a
#ベクトルを作るときはc()で指定。

#行列
mat1 <- matrix(c(1:10),nrow=2,ncol=5)
mat1 <- matrix(c(1:10),2,5)
mat2 <- matrix(c(1:10),2,byrow=T)
mat3 <- matrix(c(1,3,2,5,7,3,2,15,2),3,3,byrow=T)
#行列はmatrix()で作成。byrow=Tで要素順が変わる。

#行列の計算
mat1+mat2
mat1-mat2
mat1*mat2
mat1/mat2
mat3%*%mat3
solve(mat3)

#行列の各要素へのアクセス(参照)
p <- mat1[1,2]
q <- mat3[,2]
r <- mat1[1,c(2,5)]

#データ読み込み
df <- read.csv(“sample-data.csv”,header=T,row.names=1)

#①読み込みたいデータ名を””で囲む
#②headerは列名があるか無いか。あるときはT(TURE)、ない時はF(FALSE)
#③row.namesは行名の指定。row.names=1だと1列目が列名になる

#データのアクセス
df[,1] #1列目のデータ
df[1,] #1行目のデータ
df$年齢 #「年齢」という名の列のデータ

#データの取り出し方
#例1:男性だけを取り出す
M.dat <- subset(df,df$性別==”M”)
M.dat <- subset(df,df[,4]==”M”)
M.dat <- df[df$性==”M”,]

#例2:体重が60kg未満の人を取り出す
weight60 <- subset(df,df$体重<60)
weight60 <- df[df[,6]<60,]

#例3:男性かつ体重が60kg未満の人を取り出す
M.60 <- subset(df,df$性別==”M” & df$体重<60)
M.60 <- df[df[,4]==”M” & df[,6]<60,]

h<-subset(df,df$肺活量>=4000)
h
h1<-subset(df,df$肺活量>=3000 & df$肺活量<=4000)
h1
b<-subset(df,df$病気==1 & df$体重>=70)
b

###############
#データ読み込み
df <- read.csv(“sample-data.csv”,header=T,row.names=1)
##ディレクトリが異なる場合はファイルのアドレスも記述

#便利なコマンド
x <- numeric(nrow(df))
ncol(df)
#nrowは行数,numericは0を引数個ならべる関数(列数ならncol)

x <- rep(0,nrow(df))
#repは繰り返し

da <- 12345
da.ch <- as.character(da)
da <- as.numeric(da.ch)
is.character(da)
is.numeric(da)
#as.~は型変更、is.~は型確認(numeric,character,vector,matrix,data.frameで使われる)

y <- c(NA,1,2,NA,5)
y.na <- is.na(y)
#is.naはNAかどうかを調べる

z <- c(2,4,6,1,3,5)
length(z)
#lengthはベクトルの長さを測る

table(df$性別)
table(df$性別,df$病気)
#tableは要素のカウントをする

colnames(df)
rownames(df)
#colnames, rownamesはデータの行名、列名を見る
#ベクトルのときは行、列は関係ないのでnames()を使う

a <- c(1,2,3,1,2,3,1,2,3)
unique(a)
#uniqueは重複する要素を除いたベクトルを返す

b <- c(1,4,5,3,7,8,9,9,2,5)
mean(b)
var(b)
sd(b)
median(b)
sum(b)
max(b)
min(b)
#tips参照

rev(b)
order(b)
sort(b)
b[order(b)]
#revは順番を逆にする
#orderはベクトルの小さい順に何番目かを返す
#sortは昇順にする
#降順は?#

k <- 1:10
l <- seq(2,20,by=2)
union(k,l)
setdiff(k,l)
intersect(k,l)
#和集合、差集合、積集合

na.omit(y)
#欠損(NA)を除く

###※これらのコマンドはすぐに覚える必要はありませんが、こんなのがあるという事を頭の隅に入れておいてください(全部tipsに載っています)###
#for文とif文
ans <- 0
for(i in 1:10){ #iは1~10
if(i%%2==0){ #%%は剰余を表す
ans <- ans*i
}else{
ans <- ans-i
}
}

ans <- 0
for(j in c(2,4,6,8,10)){ #jに2,4,…というベクトルの要素を順に代入
ans <- ans+j
}
#回帰分析
y <- c(1,3,4,10,5,1,3,14,21)
x <- c(10,20,10,40,50,10,10,20,70)

ans <- lm(y~x)
s.ans <- summary(ans)
coe <- s.ans$coefficient
#coe <- coef(s.ans)やs.ans$coefも同じ出力
N <- length(y)
aic <- AIC(ans)
result <- cbind(coe,aic,N)
result[2,5:6] <- “”
write.table(result,”回帰分析.csv”,append=T,quote=F,sep=”,”,row.names=T,col.names=T)
#出力内容、ファイル名、append:ファイルに付け足しokか、quote:””を出力するか、sep:データの区切りを何にするか、row.names:行名を出力するか、col.names:列名を出力するか#

write.table(matrix(c(“”,colnames(result)),nrow=1),”回帰分析.csv”,append=T,quote=F,sep=”,”,row.names=F,col.names=F)
write.table(result,”回帰分析.csv”,append=T,quote=F,sep=”,”,row.names=T,col.names=F)
ans <- lm(df$肺活量~df$血圧)
s.ans <- summary(ans)
coe <- s.ans$coefficient
N <- nrow(df)
aic <- AIC(ans)
result <- cbind(coe,aic,N)
result[2,5:6] <- “”

write.table(matrix(c(“”,colnames(result)),nrow=1),”回帰分析.csv”,append=T,quote=F,sep=”,”,row.names=F,col.names=F)
write.table(result,”回帰分析.csv”,append=T,quote=F,sep=”,”,row.names=T,col.names=F)

################
#データ読み込み
df <- read.csv(“sample-data.csv”,header=T,row.names=1)
###重回帰分析###
ans <- lm(df$肺活量~df$血圧+df$体重)
s.ans <- summary(ans)
coe <- s.ans$coefficient
N <- nrow(df)
aic <- AIC(ans)
result <- cbind(coe,aic,N)
result[2:nrow(result),5:6] <- “”
filename=”重回帰分析.csv”
write.table(matrix(c(“”,colnames(result)),nrow=1),filename,append=T,quote=F,sep=”,”,row.names=F,col.names=F)
write.table(result,filename,append=T,quote=F,sep=”,”,row.names=T,col.names=F)
write.table(“”,filename,append=T,quote=F,sep=”,”,row.names=F,col.names=F)

dat <- df[,c(3,2,1)]
ans <- lm(df$肺活量~.,data=dat)
s.ans <- summary(ans)
coe <- s.ans$coefficient
N <- nrow(df)
aic <- AIC(ans)
result <- cbind(coe,aic,N)
result[2:nrow(result),5:6] <- “”

#excelに複数のデータの結果を付け加えるときにはappendをTにする
filename <- “重回帰出力test1.csv”
write.table(matrix(c(“”,colnames(result)),nrow=1),filename,append=T,quote=F,sep=”,”,row.names=F,col.names=F)
write.table(result,filename,append=T,quote=F,sep=”,”,row.names=T,col.names=F)
write.table(“”,filename,append=T,quote=F,sep=”,”,row.names=F,col.names=F)

#練習
#y:年齢,x1~3:血圧、性別、病気
#ダミー変数について
df$病気
factor(df$病気)#要素をグループ化

###logistic回帰分析###
dat <- df[,c(5,1,2,6)]
ans <- glm(df$病気~.,data=dat,family=binomial)
s.ans <- summary(ans)
coe <- s.ans$coefficient
N <- nrow(df)
aic <- AIC(ans)
result <- cbind(coe,aic,N)
result[2:nrow(result),5:6] <- “”

###cox回帰###
#cox回帰では応答変数としてtimeとstatusの二つを用意する。
#timeは時間、statusは死亡などのactionが1,0ではいる。
library(survival)
#データ読み込み
df <- read.table(“cox_sample_data.txt”,header=T,sep=”\t”)
dat <- df[,c(1,2,4,5,6)]
dat[,4] <- factor(dat[,4])
dat[,5] <- factor(dat[,5])
ans <- coxph(Surv(df$OS,df$死亡の有無)~.,data=dat)
s.ans <- summary(ans)
#以下ほしいデータを抽出
coef <- round(s.ans$coef,5)
aic <- round(-2*max(ans$loglik)+2*nrow(coef),5)
n <- s.ans $n

na.omit(dat)
nrow(na.omit(dat))

ci <- round(s.ans$conf.int,5)
ci.low <- ci[,3]
ci.up <- ci[,4]
result <- cbind(coef,ci.low,ci.up,aic,n)
result[2:nrow(result),8:9] <- “”
###相関係数###
a <- c(1:10)
b <- c(1,5,23,6,4,2,9,8,10,6)
c <- c(2,4,6,8,10,13,8,2,19,29)
dat <- cbind(a,b,c)

cor(a,b)
cor(dat)
###関数###
fun1 <- function(x){
a <- x^2
b <- x+10
ans <- a+b
return(ans)
}
fun.ans <- fun1(10)
df <- read.csv(“sample-data.csv”,header=T,row.names=1)
fun2 <- function(dat,num){
ans <- glm(df$病気~.,data=dat,family=binomial)
s.ans <- summary(ans)
coe <- s.ans$coefficient
N <- nrow(df)
aic <- AIC(ans)
result <- cbind(coe,aic,N)
result[2:nrow(result),5:6] <- “”

filename <- paste(“logistic回帰-sample-“,num,”変量.csv”,sep=””)
write.table(matrix(c(“”,colnames(result)),nrow=1),filename,append=F,quote=F,sep=”,”,row.names=F,col.names=F)
write.table(result,filename,append=T,quote=F,sep=”,”,row.names=T,col.names=F)
write.table(“”,filename,append=T,quote=F,sep=”,”,row.names=F,col.names=F)
}

fun2(df[,c(5,1,2,3)],3)

############
#データ読み込み
df <- read.csv(“sample-data.csv”,header=T,row.names=1)

#t検定
#対応の無い検定
a <- c(10,12,9,15,8,2)
b <- c(20,18,12,23,16)

t.ans1 <- t.test(a,b,var=T) #等分散仮定 sudent
t.ans2 <- t.test(a,b,var=F) #不等分散仮定 ウェルチ
t.ans3 <- wilcox.test(a,b) #ノンパラメトリック ウィルコクソン

#対応のある検定
c <- c(10,11,6,3,19,14)
d <- c(15,15,8,4,21,10)

t.ans1 <- t.test(c,d,paired=T)
t.ans2 <- wilcox.test(c,d,paired=T)

#30歳以上
over30 <- df[df$年齢>=30,]
#30歳未満
under30 <- df[df$年齢<30,]

#t検定をして,p値
t.fun <- function(da1,da2){
ans1 <- t.test(da1,da2,var=T)
ans2 <- t.test(da1,da2,var=F)
ans3 <- wilcox.test(da1,da2)
result <- c(ans1$p.value,ans2$p.value,ans3$p.value)
return(result)
}

ans <- c()
#各変数において平均値の差の検定
for(i in c(2,3,6)){
dat1 <- over30[,i]
dat2 <- under30[,i]
ans <- rbind(ans,t.fun(dat1,dat2))
}
rownames(ans) <- colnames(df)[c(2,3,6)]
colnames(ans) <- c(“t検定(等分散仮定)”,”t検定(不等分散仮定)”,”マンホイットニーのU検定”)
###独立性の検定###
mat1 <- matrix(c(100,70,40,20,60,90,40,10,30,20,90,60,10,10,80,110),4,4,byrow=T)
ans1 <- chisq.test(mat1)$p.value

mat2 <- matrix(c(20,40,10,50),2,2)
#カイ二乗検定、フィッシャーの正確検定
ans <- c(chisq.test(mat2)$p.value,fisher.test(mat2)$p.value)
###plot###
#色々な図のプロットがあるので詳しくはr-tipsを参照
x <- c(1,2,3,2,7,5,9,1)
y <- c(14,20,21,15,36,27,40,8)
ans <- lm(y~x)

plot(x,y,lty=6)

lines(x,fitted(ans),col=”red”)
dev.off()

plot(x,y,pch=6,col=”blue”)
lines(x,fitted(ans),lty=3,col=”red”)
#図の保存の仕方
#(直接ファイルから保存することもできるが解析量が多くなるとお勧めできない)
f1 <- “sample1.png”
f2 <- “sample2.jpeg”
png(f1, width = 800, height = 600)
plot(x,y)
lines(x,fitted(ans),col=”red”)
dev.off()

plot(x,y)
lines(x,fitted(ans),col=”red”)
dev.copy(jpeg,f2,width=900,height=675)
dev.off()

※R言語入門のトップページはこちら

(totalcount 3,880 回, dailycount 10回 , overallcount 16,405,036 回)

ライター:

R言語入門

single-banner

COMMENT

コメントを残す

メールアドレスが公開されることはありません。
*は必須項目です。




CAPTCHA