PDFにページ数を挿入する場合

PDFにページ数を挿入する場合
PDFのファイル名はNo1,No2,No3の三枚
同じフォルダに入れる

\documentclass[uplatex,11pt]{jsarticle}
\usepackage[top=15truemm,bottom=15truemm,left=15truemm,right=15truemm]{geometry}     %余白

%図
\usepackage[dvipdfmx]{graphicx}
%図の位置調整　改行後に[H]
\usepackage{here}  

% 章が進むごとに図番号をリセットする
\makeatletter% プリアンブルで定義開始
\renewcommand{\thefigure}{\arabic{figure}}
\@addtoreset{figure}{section}
\makeatother % プリアンブルで定義終了

%図とキャプションを狭くする
\setlength\abovecaptionskip{5pt}

%---------------------------------------------------%
\begin{document} 
%---------- 1 ----------%
\begin{figure}[H] %現在の位置に
  \centering
    \includegraphics[width=17cm]{No1.pdf} %図の挿入
\end{figure}
%---------- 2 ----------%
\newpage
\begin{figure}[H] %現在の位置に
  \centering
    \includegraphics[width=17cm]{No2.pdf} %図の挿入
\end{figure}
%---------- 3 ----------%
\newpage
\begin{figure}[H] %現在の位置に
  \centering
    \includegraphics[width=17cm]{No3.pdf} %図の挿入
\end{figure}

\end{document} 

層別した散布図と回帰直線

以下のようなグラフ
忘れないように記載しておきます

データサンプル

Group <- 
c(1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 
  3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 
  2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3)

pre_post <- 
structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
            2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 
            1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
            1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("治療後", "治療前"
            ), class = "factor")

Outcome <- 
c(23.5, 21.17, 15.5, 20.5, 21.17, 11.17, 19, 20.5, 20.17, 19.17, 
  11, 12, 24.83, 23.17, 22.5, 24.17, 26.5, 22, 11.83, 16.83, 10.5, 
  19.5, 19.83, 14, 12.936, 12.936, 18, 7.864, 6.264, 16.664, 20.136, 
  21.736, 3.6, 16.664, 13.2, 2.8, 12.936, 10.264, 18.8, 22.5, 9.83, 
  7.83, 20.83, 25.17, 27.17, 4.5, 20.83, 16.5, 3.5, 16.17, 12.83, 
  23.5, 21.17, 15.5)

#data に取り込みます
data <- data.frame(Group, pre_post, Outcome)

グラフを作成します

治療前 <- data[data$pre_post=="治療前",3]
治療後 <- data[data$pre_post=="治療後",3]
plot(治療前, 治療後, col=c(data$Group), pch =16, xlab = "治療前", ylab="治療後" )

#Group=1 黒のライン
data1 <- subset(data, data$Group==1)
治療前 <- data1[data1$pre_post=="治療前",3]
治療後 <- data1[data1$pre_post=="治療後",3]
xs <- seq(0, 30, length=100)
fit <- glm(治療後~治療前)
lines(xs, fit$coef[1] + xs*fit$coef[2], col=1)

#Group=2 赤のライン
data2 <- subset(data, data$Group==2)
治療前 <- data2[data2$pre_post=="治療前",3]
治療後 <- data2[data2$pre_post=="治療後",3]
xs <- seq(0, 30, length=100)
fit <- glm(治療後~治療前)
lines(xs, fit$coef[1] + xs*fit$coef[2], col=2)

#Group=3 緑のライン
data3 <- subset(data, data$Group==3)
治療前 <- data3[data3$pre_post=="治療前",3]
治療後 <- data3[data3$pre_post=="治療後",3]
xs <- seq(0, 30, length=100)
fit <- glm(治療後~治療前)
lines(xs, fit$coef[1] + xs*fit$coef[2], col=3)

ロジスティック回帰の基準変更

data2 <- data.frame(c(rep("あり",15),rep("なし",15)),
               c(rep("有効",12),rep("無効",3),rep("有効",4),"無効","有効","有効",rep("無効",8)))
colnames(data2) <- c("筋トレ","効果")
xtabs(~ 効果 + 筋トレ, data = data2)  

      筋トレ
効果   あり なし
  無効    3    9
  有効   12    6

オブジェクトの内容を情報付きで簡潔に表示 str()

str(data2)

'data.frame':  30 obs. of  2 variables:
 $ 筋トレ: Factor w/ 2 levels "あり","なし": 1 1 1 1 1 1 1 1 1 1 ...
 $ 効果  : Factor w/ 2 levels "無効","有効": 2 2 2 2 2 2 2 2 2 2 ...

# この順番で基準が決定されます
# 1番目に対する2番目の確率を求めることになります

筋トレ効果を考える場合にロジスティック回帰を実行すると、
「オッズ ＝ 有効 / 無効、オッズ比 = なし / あり」
となります

glm(効果 ~ 筋トレ, family = binomial, data=data2)

Coefficients:
(Intercept)   筋トレなし  
      1.386       -1.792  

オッズ比を確認してみます

exp(-1.792) = 0.167

#オッズ ＝ 有効 / 無効、オッズ比 = なし / あり　なので

      筋トレ
効果   あり なし
  無効    3    9
  有効   12    6

(6/9) / (12/3) = 0.167

オッズの基準を変更する場合には

data2$効果 <- relevel(data2$効果, ref="有効")
glm(効果 ~ 筋トレ, family = binomial, data=data2)

(Intercept)   筋トレなし  
     -1.386        1.792  

#オッズ ＝ 無効 / 有効  、オッズ比 = なし / あり　なので
オッズ比 = exp(1.7918) = 6.0
オッズ比 = (9/6) / (3/12) = 6.0

オッズ比の基準を変更する場合も同じです

#オッズ = 有効／無効、オッズ比 ：「筋トレあり」の確率
data2$筋トレ <- relevel(data2$筋トレ, ref="なし")
glm(効果 ~ 筋トレ, family = binomial, data=data2)

開設3周年

はてなブログさんから素敵なメッセージが届きました

表１の作り方 "tableone"

年齢   性 運動習慣
52 男 NO
55 男 NO
52 男 NO
56 男 YES
56 男 NO
49 男 NO
43 男 NO
48 男 YES
43 男 NO
42 男 YES
44 男 NO
43 男 NO
48 男 NO
51 男 NO
43 男 NO
49 男 NO
39 男 YES
38 男 NO
50 男 YES
52 女 NO
59 女 YES
65 女 NO
54 女 YES
39 女 NO
53 女 NO
39 女 NO
50 女 NO
41 女 YES
44 女 YES
56 女 NO
54 女 YES
57 女 NO
53 女 NO
51 女 NO
61 女 YES
51 女 NO
57 女 NO
59 女 YES
55 女 YES
52 女 NO

上のサンプルをコピーしてRに取り込みます

data01 <- read.table("clipboard", header = T)
install.packages("tableone")
library(tableone)

#一覧にしたい変数名を記載
vars <- c("年齢","性")

#名義変数を指定
# 注意）IDがある場合にはファクターに入れないこと
facvars <- c("性")

#strataで分割する
table <- CreateTableOne(
  vars = vars, 
  strata = "運動習慣", 
  data = data01, 
  factorVars = facvars) 

table 

p値：連続変数＝分散分析/ｔ検定 , 名義変数＝ｶｲ二乗検定

                  Stratified by 運動習慣
                   NO            YES           p      test
  n                   27            13                    
  年齢 (mean (SD)) 49.67 (6.44)  50.92 (7.48)   0.586     
  性 = 男 (%)         14 (51.9)      5 (38.5)   0.648     

性をfisherの正確検定、年齢を正規分布に従わないとしてIQRを提示する場合
性別は一行に両レベル表示 quote = TRUEでエクセルで区切り位置設定で必要

print(table, nonnormal = c("年齢"),
      exact = c("性"), cramVars = "性",
      quote = TRUE)

                     Stratified by 運動習慣
                      NO                   YES                  p      test   
  n                      27                   13                              
  年齢 (median [IQR]) 51.00 [43.50, 53.00] 54.00 [44.00, 56.00]  0.478 nonnorm
  性 = 女/男 (%)      13/14 (48.1/51.9)      8/5 (61.5/38.5)     0.511 exact  

詳細は下記を参照ください

tableone: 医学研究で必要な"表1"を作成するRパッケージ

層別カテゴリー化した変数を追加

id   age
1  35
2  55
3  53
4  56
5  66
6  49
7  43
8  48
9  64
10 42
11 44
12 74
13 48
14 61
15 43
16 49
17 39
18 38
19 55
20 43
21 48
22 48
23 49
24 49
25 43
26 46
27 49
28 46
29 47
30 43
31 43
32 48
33 50
34 41
35 51
36 45
37 50
38 47
39 44
40 68
41 45
42 45
43 49
44 75
45 55
46 77
47 42
48 57
49 48
50 58

このデータのageを4つにカテゴリー化して、列を追加します
上記のデータをコピーしてRに転送します

data <- read.table("clipboard", header = T)

package ”dplyr” にあるmutateを使用して、列を追加します

library(dplyr)
data2 <- mutate(data, age_c = 
         if_else(data$age < 50 , "50未満", 
                 if_else(data$age >= 50 & data$age < 60, "50-60", 
                         if_else(data$age >= 60 & data$age < 70, "60-70",
                                 "70以上")
                         )))

データの確認

head(data2)

  id age  age_c
1  1  35 50未満
2  2  55  50-60
3  3  78 70以上
4  4  56  50-60
5  5  66  60-70
6  6  49 50未満

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

使用するパッケージ

install.packages("ggplot2")
library(ggplot2)

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

#乱数を使用のため、実施毎にグラフは異なります
dat <- c(c(rnorm(15)),c(rnorm(15,3,1.2)),c(rnorm(15,10,2)),c(rnorm(15,17,1.5)))
pre_post <- c(rep("前", 30),rep("後", 30))
treat <- c(rep("A",15),rep("B",15),rep("A",15),rep("B",15))
da01 <- data.frame(dat, pre_post, treat);head(dataf)

            dat pre_post treat
1  0.0004020798       前     A
2  0.3236960969       前     A
3  1.0026311892       前     A
4  0.3526106439       前     A
5 -0.1307876208       前     A
6 -1.5476749466       前     A

列pre_postの配列を確認します

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

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

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

二群比較

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

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

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

(g <- ggplot(data = da02,
       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.y = "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()) #背景の設定

単回帰分析

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

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.2681     0.5873  -0.456    0.652    
pre_post後    4.0898     0.8306   4.924 3.42e-05 ***

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

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.0459     0.6852   4.445 0.000126 ***
pre_post後    1.0633     0.9690   1.097 0.281873  

単回帰の結果を書き込んでみます

(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)

重回帰分析

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

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -0.2681     0.6381  -0.420 0.676005    
pre_post後          4.0898     0.9025   4.532 3.12e-05 ***
treatB              3.3140     0.9025   3.672 0.000539 ***
pre_post後:treatB  -3.0265     1.2763  -2.371 0.021189 *  

#分散分析表
anova(fit)

Response: dat
               Df Sum Sq Mean Sq F value    Pr(>F)    
pre_post        1  99.58  99.579 16.3016 0.0001658 ***
treat           1  48.64  48.640  7.9627 0.0065965 ** 
pre_post:treat  1  34.35  34.349  5.6231 0.0211892 *  
Residuals      56 342.08   6.109  

この記事で使用したデータフレーム
コピーしてRにペースト
da01 <- read.table("clipboard", header=T)
これで再現できます

> da01
    X         dat pre_post treat
1   1 -0.29733004       前     A
2   2 -2.63812280       前     A
3   3 -0.90097072       前     A
4   4  1.06843016       前     A
5   5 -3.03846846       前     A
6   6  2.14097694       前     A
7   7  2.47494865       前     A
8   8 -0.02154341       前     A
9   9  0.56411223       前     A
10 10  2.81826063       前     A
11 11 -1.22557353       前     A
12 12 -1.96770187       前     A
13 13 -2.41952982       前     A
14 14 -2.04788520       前     A
15 15 -1.49396856       前     A
16 16  4.42978119       前     B
17 17  1.14121396       前     B
18 18  2.81408173       前     B
19 19  3.88617124       前     B
20 20  0.27836700       前     B
21 21  7.04775580       前     B
22 22  6.76956940       前     B
23 23 -1.01243943       前     B
24 24  2.28982618       前     B
25 25  2.56697055       前     B
26 26  3.92271476       前     B
27 27  0.56754099       前     B
28 28  2.55950816       前     B
29 29  5.45781910       前     B
30 30  5.97792826       前     B
31 31  3.57321119       後     A
32 32  3.84212737       後     A
33 33  1.38192000       後     A
34 34  1.89198329       後     A
35 35  5.90437940       後     A
36 36  3.16262753       後     A
37 37  1.29383105       後     A
38 38  7.08585415       後     A
39 39  1.53787689       後     A
40 40  5.82484000       後     A
41 41  4.45300261       後     A
42 42  3.94823267       後     A
43 43  2.22314781       後     A
44 44  5.21597500       後     A
45 45 -1.10466187       後     A
46 46  1.03430878       後     B
47 47  3.26798971       後     B
48 48  4.46817078       後     B
49 49  4.30615587       後     B
50 50  3.13762549       後     B
51 51  3.12041534       後     B
52 52  5.40843117       後     B
53 53  2.01223446       後     B
54 54 -2.57111800       後     B
55 55  4.01156402       後     B
56 56  6.58882849       後     B
57 57  6.94067008       後     B
58 58  3.12482569       後     B
59 59  3.81695957       後     B
60 60  8.53560520       後     B