CPR-2023-Homework

0. 準備

ライブラリを読み込む。

pacman::p_load(tidyverse, 
               foreign, 
               rms, 
               tableone, 
               ggplot2, 
               gtsummary, 
               skimr, 
               mgcv)

データを読み込む。

dat_w <- read_csv("GustoW1.csv", col_names = TRUE)
dat_4 <- read_csv("sample4.csv", col_names = TRUE)
dat_5 <- read_csv("sample5.csv", col_names = TRUE)

欠測値の有無を確認する。
一切ないことが分かる。

sum(is.na(dat_w), is.na(dat_4), is.na(dat_5))

## [1] 0

1. 課題1

新たにregionという名前の変数を作成し、gustowではwest,
sample4ではsample4, sample5ではsample5を入力してください。

まず、…1というID？の列を落とす。
また、dat_wだけage10という余計な変数があるので落とす。
さらに、課題の通りregionを追加する(west, sample4, sample5)。

dat_w <- dat_w %>% 
  select(-...1, -AGE10) %>% 
  mutate(region = "west")
dat_4 <- dat_4 %>% 
  select(-...1) %>% 
  mutate(region = "sample4")
dat_5 <- dat_5 %>% 
  select(-...1) %>% 
  mutate(region = "sample5")

# 例として表示
glimpse(dat_w)

## Rows: 2,188
## Columns: 27
## $ DAY30  <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, …
## $ SEX    <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, …
## $ AGE    <dbl> 70.313, 59.844, 59.023, 80.375, 64.750, 53.883, 48.594, 72.453,…
## $ A65    <dbl> 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ KILLIP <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, …
## $ SHO    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ DIA    <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, …
## $ HYP    <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, …
## $ HRT    <dbl> 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, …
## $ ANT    <dbl> 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, …
## $ PMI    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, …
## $ HEI    <dbl> 177.3, 172.0, 170.0, 154.9, 167.0, 171.0, 177.3, 187.0, 197.0, …
## $ WEI    <dbl> 84.0, 115.0, 76.0, 50.0, 97.4, 100.0, 82.0, 90.9, 96.4, 61.0, 8…
## $ HTN    <dbl> 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, …
## $ SMK    <dbl> 3, 1, 1, 3, 1, 1, 2, 2, 2, 2, 1, 2, 2, 1, 1, 2, 1, 3, 3, 1, 1, …
## $ LIP    <dbl> 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ PAN    <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, …
## $ FAM    <dbl> 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, …
## $ STE    <dbl> 1, 6, 3, 3, 2, 4, 3, 7, 8, 3, 5, 3, 2, 5, 4, 1, 4, 1, 3, 2, 7, …
## $ ST4    <dbl> 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, …
## $ TTR    <dbl> 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, …
## $ ESAMP  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ GRPL   <dbl> 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, …
## $ GRPS   <dbl> 5, 2, 6, 2, 6, 6, 6, 1, 1, 6, 6, 6, 6, 2, 1, 6, 6, 2, 1, 6, 4, …
## $ REGL   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ HIG    <dbl> 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, …
## $ region <chr> "west", "west", "west", "west", "west", "west", "west", "west",…

2. 課題2

gustow, sample 4, sample5をマージして、
重複がないことを確認してください。

rbindでマージする。
変数の数が変わっていないので、データセット間で異なる変数はない。

dat <- rbind(dat_w, dat_4, dat_5)

glimpse(dat)

## Rows: 3,402
## Columns: 27
## $ DAY30  <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, …
## $ SEX    <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, …
## $ AGE    <dbl> 70.313, 59.844, 59.023, 80.375, 64.750, 53.883, 48.594, 72.453,…
## $ A65    <dbl> 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ KILLIP <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, …
## $ SHO    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ DIA    <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, …
## $ HYP    <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, …
## $ HRT    <dbl> 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, …
## $ ANT    <dbl> 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, …
## $ PMI    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, …
## $ HEI    <dbl> 177.3, 172.0, 170.0, 154.9, 167.0, 171.0, 177.3, 187.0, 197.0, …
## $ WEI    <dbl> 84.0, 115.0, 76.0, 50.0, 97.4, 100.0, 82.0, 90.9, 96.4, 61.0, 8…
## $ HTN    <dbl> 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, …
## $ SMK    <dbl> 3, 1, 1, 3, 1, 1, 2, 2, 2, 2, 1, 2, 2, 1, 1, 2, 1, 3, 3, 1, 1, …
## $ LIP    <dbl> 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ PAN    <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, …
## $ FAM    <dbl> 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, …
## $ STE    <dbl> 1, 6, 3, 3, 2, 4, 3, 7, 8, 3, 5, 3, 2, 5, 4, 1, 4, 1, 3, 2, 7, …
## $ ST4    <dbl> 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, …
## $ TTR    <dbl> 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, …
## $ ESAMP  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ GRPL   <dbl> 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, …
## $ GRPS   <dbl> 5, 2, 6, 2, 6, 6, 6, 1, 1, 6, 6, 6, 6, 2, 1, 6, 6, 2, 1, 6, 4, …
## $ REGL   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ HIG    <dbl> 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, …
## $ region <chr> "west", "west", "west", "west", "west", "west", "west", "west",…

3. 課題3

連続変数はどれでしょうか？
連続変数は、ヒストグラム、箱ひげ図、などを描いて、分布を確認してください。
注：以後の解析では、連続変数以外は、factor変数として、Rに指定して解析を実施するようにしてください。

glimpse(dat)から、AGE, HEI, WEIが連続変数であることが分かる。
それ以外をfactor型に指定する。

dat <- dat %>% 
  mutate(across(.cols = c(-AGE, -HEI, -WEI), .fns = factor))

AGEのヒストグラム

hist_AGE <- ggplot(data = dat) + 
  geom_histogram(aes(x = AGE), bins = 30) + 
  theme_classic()

hist_AGE

DAY30でカテゴリ分けしたAGEの箱ひげ図

box_AGE <- ggplot(data = dat) + 
  geom_boxplot(aes(x = DAY30, y = AGE)) + 
  theme_classic()

box_AGE

HEIについても同様に（見やすさを考えてbins=20）

hist_HEI <- ggplot(data = dat) + 
  geom_histogram(aes(x = HEI), bins = 20) + 
  theme_classic()

hist_HEI

box_HEI <- ggplot(data = dat) + 
  geom_boxplot(aes(x = DAY30, y = HEI)) + 
  theme_classic()

box_HEI

WEIについても同様に（bins=20）

hist_WEI <- ggplot(data = dat) + 
  geom_histogram(aes(x = WEI), bins = 20) + 
  theme_classic()

hist_WEI

box_WEI <- ggplot(data = dat) + 
  geom_boxplot(aes(x = DAY30, y = WEI)) + 
  theme_classic()

box_WEI

4. 課題4

論文のTable 1のようにOverall populationとDeathsの集団の列を作成し、
基本属性を表す表を作成してください。

gtsummaryのtbl_summaryを使用してみる。

tbl_kadai4 <- dat %>% 
  tbl_summary(by = DAY30, 
              percent = "row") %>% 
  add_n() %>% 
  add_overall() %>% 
  modify_header(label ~ "**Variable**") %>%
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Death Before DAY30**") %>% 
  modify_caption("**Kadai4. Baseline Characteristics**") %>%
  bold_labels()

tbl_kadai4

**Kadai4. Baseline Characteristics**
Variable	N	Overall, N = 3,402¹	Death Before DAY30
Variable	N	Overall, N = 3,402¹	0, N = 3,191¹	1, N = 211¹
SEX	3,402
0		2,541 (100%)	2,413 (95%)	128 (5.0%)
1		861 (100%)	778 (90%)	83 (9.6%)
AGE	3,402	61 (51, 70)	60 (51, 69)	73 (65, 78)
A65	3,402
0		2,079 (100%)	2,025 (97%)	54 (2.6%)
1		1,323 (100%)	1,166 (88%)	157 (12%)
KILLIP	3,402
1		2,927 (100%)	2,798 (96%)	129 (4.4%)
2		415 (100%)	361 (87%)	54 (13%)
3		42 (100%)	24 (57%)	18 (43%)
4		18 (100%)	8 (44%)	10 (56%)
SHO	3,402
0		3,342 (100%)	3,159 (95%)	183 (5.5%)
1		60 (100%)	32 (53%)	28 (47%)
DIA	3,402
0		2,949 (100%)	2,776 (94%)	173 (5.9%)
1		453 (100%)	415 (92%)	38 (8.4%)
HYP	3,402
0		3,104 (100%)	2,935 (95%)	169 (5.4%)
1		298 (100%)	256 (86%)	42 (14%)
HRT	3,402
0		2,312 (100%)	2,209 (96%)	103 (4.5%)
1		1,090 (100%)	982 (90%)	108 (9.9%)
ANT	3,402
0		2,154 (100%)	2,064 (96%)	90 (4.2%)
1		1,248 (100%)	1,127 (90%)	121 (9.7%)
PMI	3,402
0		2,822 (100%)	2,674 (95%)	148 (5.2%)
1		580 (100%)	517 (89%)	63 (11%)
HEI	3,402	173 (165, 178)	173 (165, 178)	170 (160, 175)
WEI	3,402	80 (70, 90)	80 (70, 91)	71 (61, 83)
HTN	3,402
0		2,056 (100%)	1,943 (95%)	113 (5.5%)
1		1,346 (100%)	1,248 (93%)	98 (7.3%)
SMK	3,402
1		1,349 (100%)	1,297 (96%)	52 (3.9%)
2		1,107 (100%)	1,027 (93%)	80 (7.2%)
3		946 (100%)	867 (92%)	79 (8.4%)
LIP	3,402
0		2,063 (100%)	1,929 (94%)	134 (6.5%)
1		1,339 (100%)	1,262 (94%)	77 (5.8%)
PAN	3,402
0		2,228 (100%)	2,114 (95%)	114 (5.1%)
1		1,174 (100%)	1,077 (92%)	97 (8.3%)
FAM	3,402
0		1,831 (100%)	1,704 (93%)	127 (6.9%)
1		1,571 (100%)	1,487 (95%)	84 (5.3%)
STE	3,402
0		52 (100%)	49 (94%)	3 (5.8%)
1		133 (100%)	126 (95%)	7 (5.3%)
2		403 (100%)	385 (96%)	18 (4.5%)
3		1,066 (100%)	1,026 (96%)	40 (3.8%)
4		490 (100%)	458 (93%)	32 (6.5%)
5		442 (100%)	407 (92%)	35 (7.9%)
6		373 (100%)	343 (92%)	30 (8.0%)
7		297 (100%)	265 (89%)	32 (11%)
8		122 (100%)	112 (92%)	10 (8.2%)
9		18 (100%)	15 (83%)	3 (17%)
10		4 (100%)	3 (75%)	1 (25%)
11		2 (100%)	2 (100%)	0 (0%)
ST4	3,402
0		2,144 (100%)	2,044 (95%)	100 (4.7%)
1		1,258 (100%)	1,147 (91%)	111 (8.8%)
TTR	3,402
0		1,412 (100%)	1,346 (95%)	66 (4.7%)
1		1,990 (100%)	1,845 (93%)	145 (7.3%)
ESAMP	3,402
1		3,402 (100%)	3,191 (94%)	211 (6.2%)
GRPL	3,402
1		1,083 (100%)	1,012 (93%)	71 (6.6%)
2		1,534 (100%)	1,446 (94%)	88 (5.7%)
4		785 (100%)	733 (93%)	52 (6.6%)
GRPS	3,402
1		400 (100%)	379 (95%)	21 (5.2%)
2		259 (100%)	239 (92%)	20 (7.7%)
3		282 (100%)	260 (92%)	22 (7.8%)
4		329 (100%)	304 (92%)	25 (7.6%)
5		858 (100%)	810 (94%)	48 (5.6%)
6		489 (100%)	466 (95%)	23 (4.7%)
9		327 (100%)	313 (96%)	14 (4.3%)
10		195 (100%)	175 (90%)	20 (10%)
11		263 (100%)	245 (93%)	18 (6.8%)
REGL	3,402
1		2,617 (100%)	2,458 (94%)	159 (6.1%)
2		785 (100%)	733 (93%)	52 (6.6%)
HIG	3,402
0		1,771 (100%)	1,710 (97%)	61 (3.4%)
1		1,631 (100%)	1,481 (91%)	150 (9.2%)
region	3,402
sample4		785 (100%)	733 (93%)	52 (6.6%)
sample5		429 (100%)	405 (94%)	24 (5.6%)
west		2,188 (100%)	2,053 (94%)	135 (6.2%)
¹ n (%); Median (IQR)

5. 課題5

各変数（‘AGE’,‘SEX’,‘DIA’,‘HYP’,‘HRT’,‘ANT’,‘PMI’,‘HTN’,‘SMK’,‘LIP’,‘PAN’,‘FAM’,‘TTR’,‘STE’,‘WEI’,‘HEI’）とDAY30との間の2変量の関連を推定し、カイ2乗統計量、p値、自由度、オッズ比、95%信頼区間を示してください。

結果を投入する空のリストを作成する。

logreg_result <- list()
anova_result <- list()

推定値をまとめる空のベクトルを作成する。

ChiStat <- Pvalue <- Df <- OR <- CIlow <- CIhigh <- c()

DAY30との関連を調べる変数をまとめる。

vars <- c("AGE", "SEX", "DIA", "HYP", "HRT", "ANT", "PMI", "HTN", "SMK", "LIP", "PAN", "FAM", "TTR", "STE", "WEI", "HEI")

変数ごとのコード実行回数が多いので、forでまとめる。
コメントアウトにstep-by-stepの手順を記載した。

for (i in vars){
  # varsそれぞれについて、DAY30を目的変数としたロジスティック回帰分析を実行する
  logreg_formula <- as.formula(paste("DAY30 ~ ", i))
  logreg_model <- glm(logreg_formula, data = dat, family = binomial)
  
  # それぞれのモデルをリストに格納する。リスト内の要素名＝変数名（AGE, SEX, ...）とした
  logreg_result[[i]] <- logreg_model
  
  # さらに、別のリストにanovaの結果を格納する
  anova_result[[i]] <- anova(logreg_model, test = "Chisq")
  
  # 興味のある変数を取り出す（カイ二乗統計量、p値、...）
  ChiStat <- c(ChiStat, anova_result[[i]][["Deviance"]][2])
  Pvalue <- c(Pvalue, anova_result[[i]][["Pr(>Chi)"]][2])
  Df <- c(Df, anova_result[[i]][["Df"]][2])
  OR <- c(OR, exp(logreg_result[[i]][["coefficients"]][2]))
  CIlow <- c(CIlow, exp(confint(logreg_result[[i]])[2, 1]))
  CIhigh <- c(CIhigh, exp(confint(logreg_result[[i]])[2, 2]))
}

# 結果を一つの表にまとめる
result <- tibble(Variables = vars, 
                 ChiStat = ChiStat, 
                 Pvalue = Pvalue, 
                 Df = Df, 
                 OR = OR, 
                 CIlow = CIlow, 
                 CIhigh = CIhigh)

result

## # A tibble: 16 × 7
##    Variables ChiStat   Pvalue    Df    OR CIlow CIhigh
##    <chr>       <dbl>    <dbl> <int> <dbl> <dbl>  <dbl>
##  1 AGE       191.    2.01e-43     1 1.10  1.08   1.12 
##  2 SEX        21.4   3.69e- 6     1 2.01  1.50   2.68 
##  3 DIA         3.96  4.67e- 2     1 1.47  1.01   2.10 
##  4 HYP        27.1   1.89e- 7     1 2.85  1.96   4.05 
##  5 HRT        35.4   2.70e- 9     1 2.36  1.78   3.12 
##  6 ANT        39.6   3.13e-10     1 2.46  1.86   3.27 
##  7 PMI        22.6   2.00e- 6     1 2.20  1.61   2.99 
##  8 HTN         4.39  3.62e- 2     1 1.35  1.02   1.79 
##  9 SMK        23.4   8.39e- 6     2 1.94  1.36   2.79 
## 10 LIP         0.781 3.77e- 1     1 0.878 0.655  1.17 
## 11 PAN        12.6   3.89e- 4     1 1.67  1.26   2.21 
## 12 FAM         3.70  5.43e- 2     1 0.758 0.569  1.01 
## 13 TTR         9.98  1.58e- 3     1 1.60  1.19   2.17 
## 14 STE        33.1   5.08e- 4    11 0.907 0.242  4.34 
## 15 WEI        48.4   3.40e-12     1 0.969 0.960  0.978
## 16 HEI        31.1   2.49e- 8     1 0.962 0.949  0.975

6. 課題6

目的変数をDEATH、予測変数を’AGE’,‘SEX’,‘DIA’,‘HYP’,‘HRT’,‘ANT’,‘PMI’,‘HTN’,‘SMK’,‘LIP’,‘PAN’,‘FAM’,‘TTR’,‘STE’,‘WEI’,’HEI’として、ロジスティック回帰モデルに投入することを考えているとします。
全ての変数を投入したときに多重共線性の問題がなさそうか確認してください。

モデルを作成する。

model_6_1 <- glm(DAY30 ~ AGE + SEX + DIA + HYP + HRT + ANT + PMI + HTN + SMK + LIP + PAN + FAM + TTR + STE + WEI + HEI, 
                 data = dat, 
                 family = binomial)

vifを確認する→STEのカテゴリで一部10を超えているものがある。

vif(model_6_1)

##       AGE      SEX1      DIA1      HYP1      HRT1      ANT1      PMI1      HTN1 
##  1.403436  2.367537  1.093940  1.089308  1.089584  1.441947  1.201310  1.118969 
##      SMK2      SMK3      LIP1      PAN1      FAM1      TTR1      STE1      STE2 
##  1.749988  1.810190  1.099194  1.207484  1.081687  1.018052  3.552507  7.189617 
##      STE3      STE4      STE5      STE6      STE7      STE8      STE9     STE10 
## 13.190759 11.022842 11.494863 10.367469 10.626031  4.514602  1.794085  1.307022 
##     STE11       WEI       HEI 
##  1.000004  1.809314  2.791547

STEなしのモデルを再度作成する。

model_6_2 <- glm(DAY30 ~ AGE + SEX + DIA + HYP + HRT + ANT + PMI + HTN + SMK + LIP + PAN + FAM + TTR + WEI + HEI, 
                 data = dat, 
                 family = binomial)

vifを確認する。すべて5以下であり、よさそう。

vif(model_6_2)

##      AGE     SEX1     DIA1     HYP1     HRT1     ANT1     PMI1     HTN1 
## 1.371435 2.340487 1.078732 1.079873 1.077450 1.045379 1.174487 1.102458 
##     SMK2     SMK3     LIP1     PAN1     FAM1     TTR1      WEI      HEI 
## 1.747974 1.786965 1.092686 1.169744 1.069086 1.010015 1.795212 2.764134

7. 課題7

相関が強く、モデルに同時に投入すべきでないと考えられる変数はありますか？

上記の通り、STEを除いた方がよいと考える。

8. 課題8

連続変数の’AGE’,‘WEI’,’HEI’は直線的に死亡のlog(odds)と関連していますか？それとも、非直線的に関連していますか？
Restrictic cubic spline曲線を描き、確認してください。

datadist関数でデータセットの変数の範囲を定義

ddist <- datadist(dat)
options(datadist = "ddist")

課題7の結果に基づき、STEなしのデータセットで以後検討する。

dat_dat <- dat %>% 
  select(-STE)

AGE, WEI, HEIそれぞれknot 4でRCS曲線を描く。
結果からは、WEIは非直線的な関係がありそう。

fit_AGE <- lrm(DAY30 ~ rcs(AGE, 4), data = dat_dat)
plot(Predict(fit_AGE, AGE))

fit_WEI <- lrm(DAY30 ~ rcs(WEI, 4), data = dat_dat)
plot(Predict(fit_WEI, WEI))

fit_HEI <- lrm(DAY30 ~ rcs(HEI, 4), data = dat_dat)
plot(Predict(fit_HEI, HEI))

9. 課題9

’AGE’と’SEX’は’DAY30’との関連において交互作用がありそうですか。
尤度比検定から判断してください。

lrmで実施してみる。

fit_noint <- lrm(DAY30 ~ AGE + SEX, data = dat_dat)
fit_int <- lrm(DAY30 ~ AGE + SEX + AGE:SEX, data = dat_dat)

lrtest(fit_noint, fit_int)

## 
## Model 1: DAY30 ~ AGE + SEX
## Model 2: DAY30 ~ AGE + SEX + AGE:SEX
## 
## L.R. Chisq       d.f.          P 
##  1.0128974  1.0000000  0.3142097

glmで実施してみる。

fit_noint_2 <- glm(DAY30 ~ AGE + SEX, data = dat_dat, family = binomial)
fit_int_2 <- glm(DAY30 ~ AGE + SEX + AGE:SEX, data = dat_dat, family = binomial)
anova(fit_noint_2, fit_int_2, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: DAY30 ~ AGE + SEX
## Model 2: DAY30 ~ AGE + SEX + AGE:SEX
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1      3399     1389.3                     
## 2      3398     1388.3  1   1.0129   0.3142

いずれの結果からも、交互作用はなさそう。
交互作用項を入れることで、逸脱度は有意に変わるとはいえない。

CPR-2023-Homework_1

Hiroyuki Kuroda

2023-05-07