理学療法士がまとめたノート

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

単回帰分析を一括で実行する方法

tokei.net 全人類が分かる統計学より
to-kei.net

上記サイトに自分なりの注釈をつけました。
sampleもまったく同じです。

番号 <- c(1:30)
年齢 <- c(22,23,24,25,27,28,28,29,30,31,32,32,33,33,34,36,37,37,38,39,40,42,46,49,50,53,56,58,64,65)
血圧 <- c(110,128,104,112,108,126,126,104,125,120,116,124,106,134,128,128,116,132,134,116,120,130,126,140,156,124,118,144,142,144)
肺活量 <- c(110,128,104,112,108,126,126,104,125,120,116,124,106,134,128,128,116,132,134,116,120,130,126,140,156,124,118,144,142,144)
性別 <- c("M","M","F","F","M","F","F","F","F","F","M","M","F","F","M","M","M","M","F","M","F","F","M","F","M","M","M","M","F","F")
病気 <- c(1,1,0,0,0,0,1,1,1,1,1,0,0,0,1,1,1,1,0,1,1,1,0,0,1,1,1,0,1,0)
体重 <- c(79,65,53,45,80,50,43,55,47,49,64,61,48,41,70,55,70,90,39,86,50,65,58,45,60,71,62,51,40,42)
data <- data.frame(番号, 年齢, 血圧, 肺活量, 性別, 病気, 体重)

head(data, n=3)

dat <- data
names(dat)[6] <- "Y" # 6列目をYに変更

#まず関数を作成して、どんな解析の、どこの結果を記載するのか決めます
#func01は、関数名なのでなんでもOK、他は特に変更する必要なし

func01 <- function(dat,filename,response) 
{
  ans <- lm(dat$Y~.,data=dat)    #重回帰
  s.ans=summary(ans)     #重回帰の要約
  coe <- s.ans$coef    #重回帰の結果(係数、標準誤差、t値、p値)
  CI.low <- coe[,1]-1.96*coe[,2]    #係数の検定結果 95%信頼区間下限
  CI.up <- coe[,1]+1.96*coe[,2]    #係数の検定結果 95%信頼区間上限
  aic <- AIC(ans)     #AIC
  N <- nrow(dat)    #セットの行数
  result <- cbind(coe,CI.low,CI.up,aic,N)    #95%信頼区間とAICを合体
  result[2:nrow(result),c(7,8)]=""   #2行目以降の7、8列目を空白にする・・・?
  res <- paste("Y = ",response)    #単回帰の結果の1行目の文字列になる
  #ファイルの構成
  write.table(res,filename,append=T,quote=F,sep=",",row.names=F,col.names=F)
  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)
}

filename <- "てすとファイル1.csv"    #ファイル名を作成

df1 <- dat[,-1]    #番号をカット(入れてても問題なし)
for(i in 1:ncol(df1))    #列数分繰り返す(6回)
{
  if(i != 5) #5回目(5行目)を回避
  {
    #単回帰を繰り返すために2列のデータフレームを作成(ここでは5列目が目的変数)
    dat2 <- df1[ , c(5,i)]
    colnames(dat2)[1] = "Y"      #2列になったdatの1列目が目的変数
    func01(dat2, filename, "病気")    #各解析結果の1行目Y =の後に目的変数名「病気」を挿入
  }
}

#if(i != 5)
#5回目(5行目)を回避しなくてもやれる!

dat <- data
names(dat)[6] <- "Y"     # 6列目をYに変更
filename <- "てすとファイル2.csv"
df2 <- dat[,-1]     #IDをカット 入れてても大丈夫

for(i in 1:ncol(df2))    #列数分繰り返す(6回)
{
  dat3 <- df2[,c(5,i)]
  colnames(dat3)[1] = "Y" #datが2列
  func01(dat3, filename, "病気") #ここで目的変数を挿入
}

#同じ変数同士で、fitし過ぎているので注意されます
#essentially perfect fit!!!

エラーバー付きのグラフ ggplot2

更新日2020.3.13

使用するパッケージ

install.packages("ggplot2")
install.packages("ggsignif")

データセットを作成します

dat <- c(-0.29733004,-2.63812280,-0.90097072,1.06843016,-3.03846846,2.14097694,2.47494865,-0.02154341,0.56411223,2.81826063,-1.22557353,-1.96770187,-2.41952982,-2.04788520,
-1.49396856,4.42978119,1.14121396,2.81408173,3.88617124,0.27836700,7.04775580,6.76956940,-1.01243943,2.28982618,2.56697055,3.92271476,0.56754099,2.55950816,5.45781910,
5.97792826,3.57321119,3.84212737,1.38192000,1.89198329,5.90437940,3.16262753,1.29383105,7.08585415,1.53787689,5.82484000,4.45300261,3.94823267,2.22314781,5.21597500,
-1.10466187,1.03430878,3.26798971,4.46817078,4.30615587,3.13762549,3.12041534,5.40843117,2.01223446,-2.57111800,4.01156402,6.58882849,6.94067008,3.12482569,3.81695957,
8.53560520)
pre_post <- c(rep("前", 30),rep("後", 30))
treat <- c(rep("A",15),rep("B",15),rep("A",15),rep("B",15))
dataf <- data.frame(dat, pre_post, treat);head(dataf)

            dat pre_post treat
1 -0.2973300       前     A
2 -2.6381228       前     A
3 -0.9009707       前     A
4  1.0684302       前     A
5 -3.0384685       前     A
6  2.1409769       前     A

列pre_postの配列を確認します

levels(dataf$pre_post)
[1] "後" "前"

このままではx軸が後→前の配列になるので入れ替えます

da01 <- transform(dataf, pre_post= factor(pre_post, levels = c("前", "後")))
levels(da01$pre_post)
[1] "前" "後"

二群比較
パッケージggplot2を使用

library(ggplot2)
ggplot(da01, aes(x = pre_post, y = dat)) + geom_point()+theme_bw()

f:id:yoshida931:20200313164433p:plain:w500

標準誤差をエラーバーとしたグラフを描きます

#position_jitterdodgeで図の重なり具合を調整しておきます
pos_1 <- position_jitterdodge(
  jitter.width  = 0,#横にずらす
  jitter.height = 0,#縦にずらす
  dodge.width   = 0.5#カテゴリで分ける
)

(g <- ggplot(data = da01,
       aes(
         x      = pre_post,
         y      = dat,
         colour = treat,
       )) +
  geom_jitter(alpha = 0.4, position = pos_1) +  #alpha重なりの濃さ
  stat_summary(aes(x = as.numeric(pre_post) + 0.07), fun = "mean", geom = "point", size = 3, position = pos_1)+#平均値のプロット
  stat_summary(aes(x = as.numeric(pre_post) + 0.07),fun.data = "mean_se",#mean_seで標準誤差、mean_cl_normalで95%信頼区間
                  # 標準偏差   fun.ymin = function(x) mean(x) - sd(x),fun.ymax = function(x) mean(x) + sd(x),
               geom = "errorbar",
               width = 0.1,#ひげの横幅
               lwd = 1,#ひげの棒の幅
               #fun.args = list(conf.int = 0.95),#信頼区間のときに使用する
               position = pos_1
  ) +
  theme_bw()) #背景の設定

f:id:yoshida931:20200313164648p:plain:w500

単回帰分析

fitA <- lm(dat ~ pre_post , data =da01[da01$treat=="A",])
summary(fitA)

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.4656     0.5305  -0.878    0.388    
pre_post後    3.8146     0.7503   5.084  2.2e-05 ***

fitB <- lm(dat ~ pre_post , data =da01[da01$treat=="B",])
summary(fitB) 

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.2465     0.6520   4.979 2.94e-05 ***
pre_post後    0.5671     0.9221   0.615    0.544    

単回帰の結果を書き込んでみます
パッケージ ggsignifを使用

library("ggsignif")
(g2 <- g + 
geom_signif(stat = "identity",
            data = data.frame(x = c(0.87, 1.12),
                              xend = c(1.87, 2.12),
                              y = c(9, 12),
                              annotation = c("**", "N.S.")),
            aes(x = x, xend = xend, y = y, yend = y, annotation = annotation),
            col = c("red","#20B2AA")) ) +
  ylim (-5,13)

f:id:yoshida931:20200313165330p:plain:w500

重回帰分析

fit <- lm(dat ~ pre_post + treat + pre_post * treat, data =da01)
summary(fit)

                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -0.4656     0.5944  -0.783  0.43671    
pre_post後          3.8146     0.8406   4.538 3.05e-05 ***
treatB              3.7121     0.8406   4.416 4.64e-05 ***
pre_post後:treatB  -3.2475     1.1888  -2.732  0.00841 ** 

#分散分析表
anova(fit)

               Df  Sum Sq Mean Sq F value    Pr(>F)    
pre_post        1  71.995  71.995 13.5855 0.0005159 ***
treat           1  65.416  65.416 12.3440 0.0008831 ***
pre_post:treat  1  39.549  39.549  7.4629 0.0084092 ** 
Residuals      56 296.767   5.299   

簡単なエラーバーの描き方

たぶんこれが最も簡単な書き方だと思います

Aの平均 = 2
Bの平均 = 3
Aの標準偏差 = 0.4  #ここに標本標準偏差や不偏標準偏差を代入します
Bの標準偏差 = 0.7
dplot <- plot(c(1,2), c(Aの平均, Bの平均), ylim = c(0,5), xlim = c(0.5, 2.5))
arrows(1:2, c(Aの平均, Bの平均) - c(Aの標準偏差, Bの標準偏差), 1:2,  c(Aの平均, Bの平均) + c(Aの標準偏差, Bの標準偏差), lwd=1.0,angle=90,length=0.1,code=3)

f:id:yoshida931:20200313154424p:plain:w400

データセットから描く場合は、tapply関数で一括処理します

id <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
治療 <- c("A","B","A","B","B","A","A","A","B","B","A","B","A","A","A","B","A","B","B","B")
治療前BP <- c(159,168,140,170,168,135,145,150,158,170,141,159,145,141,149,165,160,168,178,175)
治療後BP <- c(145,160,135,150,157,132,140,139,150,155,130,148,138,140,142,155,145,150,160,158)
BP変化量 <- c(治療後BP - 治療前BP)

meanAB <- tapply(BP変化量,治療,mean)
sdAB <- tapply(BP変化量, 治療, sd)   #不偏標準偏差
dplot <- plot(1:2, meanAB, ylim = c(-25,10), xlim = c(0.5, 2.5))
arrows(1:2, meanAB - sdAB, 1:2, meanAB + sdAB, lwd=1.0,angle=90,length=0.1,code=3)

f:id:yoshida931:20200313155819p:plain:w400

有意差検定の結果

x1 <- 1 # x=Aの位置
x2 <- 2 # x=Bの位置
y1 <- 4 # yの位置
yhA <- y1 - ( meanAB[1] + sdAB[1] ) - 3 # x=A,y=y1から下に伸ばす距離(上限から2メモリ離す)
yhB <- y1 - ( meanAB[2] + sdAB[2] ) - 3
xx <- segments(x0=x1, y0=y1, x1=x2) 
text((x1+x2)/2, 5, paste("*"), cex=2)  #アスタリスクを書く位置
segments(x1, y1, x1, y1-yhA)
segments(x2, y1, x2, y1-yhB)

f:id:yoshida931:20200313155910p:plain:w400

比較した箱ひげにアスタリスクを入れる

アスタリスク 比較したラインを箱ひげの上限から2メモリ離した位置まで伸ばす・・・

治療 <- c("A","A","A","A","A","A","A","A","A","A","B","B","B","B","B","B","B","B","B","B")
治療前BP <- c(140,135,145,141,145,159,141,150,149,160,170,158,165,170,168,159,168,168,178,175)
治療後BP <- c(135,132,140,130,138,145,140,139,142,145,155,150,155,150,160,148,150,157,160,158)
BP変化量 <- c(治療後BP - 治療前BP)

bp <- boxplot(data$BP変化量~data$治療, ylim=c(-25,5))
#アスタリスクの下のライン
x1 <- 1 # x=Aの位置
x2 <- 2 # x=Bの位置
y1 <- 2 # yの位置
yh <- y1 - bp$stats[5,1] - 2 # x=A,y=y1から下に伸ばす距離(上限から2メモリ離す)
xx <- segments(x0=x1, y0=y1, x1=x2) 
text(1.5, 3, paste("**"),cex = 1.5)  #アスタリスクを書く位置
segments(x1, y1, x1, y1-yh)
segments(x2, y1, x2, y1-(yh + bp$stats[5,1]-bp$stats[5,2]))

[:plai]

sample

解析の練習で使うデータサンプル集 随時更新

sample1

症例 <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
性別 <- c("男性","男性","女性","男性","男性","女性","男性","男性","男性","男性","女性","女性","女性","女性","男性","女性","女性","女性","女性","男性")
治療 <- c("B","A","B","A","A","B","B","A","A","B","A","B","B","B","A","B","B","B","A","A")
治療前BP <- c(160,135,177,141,142,155,175,145,149,155,135,156,170,150,138,160,150,180,155,140)
治療後BP <- c(148,133,160,138,139,145,157,142,142,145,132,145,157,142,134,150,140,160,149,136)
BP変化量 <- c(12,2,17,3,3,10,18,3,7,10,3,11,13,8,4,10,10,20,6,4)
sample1 <- data.frame(症例,性別,治療,治療前BP,治療後BP,BP変化量 )

sample2
対応のある2群

ID <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
pre <- c(18, 8, 14, 15, 16, 13, 23, 20, 14, 20)
post <- c(23, 27, 15, 22, 15, 20, 20, 20, 18, 25)
sample2 <- data.frame(ID, pre, post)

sample3
独立した2群

グループ <- c("A","A","B","B","A","B","A","B","B","A","A","A","B","A","B","B","B","B","A","A")
ROMT <- c(18, 8, 23, 27, 14, 15, 15, 22, 15, 16, 13, 23, 20, 20, 20, 20, 18, 25, 14, 20)
sample3 <- data.frame(グループ, ROMT)

sample4
クロス表

sample4
判定 <- c(rep("非該当",30), rep("要支援", 25))
歩行 <- c(rep("可能", 20), rep("不可能", 10), rep("可能", 10), rep("不可能", 15))
sample4 <- data.frame(判定, 歩行)

テーブルをベクトルに変換

   症例 性別 治療 治療前BP 治療後BP BP変化量
1     1 男性    B      160      148       12
2     2 男性    A      135      133        2
3     3 女性    B      177      160       17
4     4 男性    A      141      138        3
5     5 男性    A      142      139        3
6     6 女性    B      155      145       10
7     7 男性    B      175      157       18
8     8 男性    A      145      142        3
9     9 男性    A      149      142        7
10   10 男性    B      155      145       10
11   11 女性    A      135      132        3
12   12 女性    B      156      145       11
13   13 女性    B      170      157       13
14   14 女性    B      150      142        8
15   15 男性    A      138      134        4
16   16 女性    B      160      150       10
17   17 女性    B      150      140       10
18   18 女性    B      180      160       20
19   19 女性    A      155      149        6
20   20 男性    A      140      136        4

上記テーブルの性別を切り出し
テキスト形式のベクトルに変換

data <- read.table("clipboard", header = T)
性別2 <- factor(data$性別)

install.packages("stringr")
library(stringr)

性別3 <- str_replace_all(性別2, pattern="k", replacement="*")
#性別3をwordで置換 スペースを,に変換

性別 <-c("男性","男性","女性","男性","男性","女性","男性","男性","男性","男性",
       "女性","女性","女性","女性","男性","女性","女性","女性","女性","男性")