BMI, hypertension in Palau, restricted cubic spline method

解析方法

この論文を参考にして，メソッドの部分を書いた方がよろしいと思います．

Restricted cubic splines were used to detect the possible nonlinear dependency of the relationship between hypertension and BMI levels, using 4 knots at prespecified locations according to the precentiles of the distribution of BMI, the 5th, 25th, 75th, and 95% percentiles.¹ The aforementioned dose-response analyses were carried out using R with package of “rms”.

グラフの下に，caption をこう書きましょう：

Association between BMI and hypertension in men/women, allowing for nonlinear effects, with 95% CIs. Model fitted with 4 knots restricted cubic spline for BMI adjusting for age, smoking, alcohol drinking, betel nut chewing, and fruit and vegetable intake frequency. Curves show ORs compared with the chosen reference BMI of 23 kg/m². BMI: Body mass index; CI, confidence interval; OR, odds ratio.

library(rms)
# 最も重要なパッケージであり，パッケージの引用方法は：
citation("rms")

To cite package ‘rms’ in publications use:

  Frank E Harrell Jr (2017). rms: Regression Modeling Strategies. R package
  version 5.1-0. https://CRAN.R-project.org/package=rms

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {rms: Regression Modeling Strategies},
    author = {Frank E {Harrell Jr}},
    year = {2017},
    note = {R package version 5.1-0},
    url = {https://CRAN.R-project.org/package=rms},
  }

ATTENTION: This citation information has been auto-generated from the package
DESCRIPTION file and may need manual editing, see ‘help("citation")’.

library(epicalc)
use(BP1)
BP1$BMI <- relevel(BP1$BMI, "23 ~ 24.9")
BP1$Smoking <- as.factor(BP1$Smoking)
BP1$Ever_drinker[is.na(BP1$Ever_drinker)] <- "unknown"
BP1$Ever_drinker <- as.factor(BP1$Ever_drinker)
BP1$Fruit <- as.factor(BP1$Fruit)
BP1$Veg <- as.factor(BP1$Veg)
BP_men <- subset(BP1, Sex == "Men")
BP_women <- subset(BP1, Sex  == "Women")
# Men
ddist <- datadist(BP_men)
options(datadist='ddist')
k <- with(BP_men, quantile(Body.Mass.Index, c(.05, 0.25, .75, .95)))
k 　#　男性で使っていたknot の位置，とその値．

      5%      25%      75%      95% 
20.12986 24.34563 32.25225 40.24364 

idx_model <- lrm(Hyt2g_medi_inclu ~ rcs(Body.Mass.Index, k)+Age+Smoking+Ever_drinker+
                   Current_Betel_Chewing + Fruit + Veg,
                 data=BP_men)
#ここは単純にロジスティクス回帰モデルにより当てはまる．
plot(Predict(idx_model, Body.Mass.Index))

ddist$limits["Adjust to","Body.Mass.Index"] <- 23
#BMI=23に対する比較するため，reference pointの設定．
idx_model <- update(idx_model)
#単純グラフの場合．
plot(Predict(idx_model, Body.Mass.Index,ref.zero=TRUE, fun=exp))

dataplot <- Predict(idx_model,Body.Mass.Index, ref.zero = TRUE, fun=exp)
#ここからは加工したグラフのプログラム：
ggplot(dataplot,aes(Body.Mass.Index, yhat)) +
#  geom_line(colour="Black", linetype="dashed", size=1.5)+
  theme(plot.subtitle = element_text(vjust = 1),
    plot.caption = element_text(vjust = 1),
    axis.line = element_line(size = 0.5,
        linetype = "solid"), panel.grid.major = element_line(colour = "gray98"),
    axis.title = element_text(size = 15),
    panel.background = element_rect(fill = "gray99",
        colour = "white", linetype = "twodash"),
    plot.background = element_rect(fill = "white")) +
  scale_x_continuous(breaks = seq(15,50,by=2.5),expression(paste("Body Mass Index", ", ", kg/m^{2})))+
  scale_y_continuous(limits = c(0,6), breaks = seq(0,6,by=0.5),"Odds Ratios (95%CI)")+
  labs(caption = NULL)+
  annotate("text", x=20, y=5.5, parse = TRUE,
           label="Men",
           size=5)+
  geom_hline(yintercept =1, linetype="dashed")

# Women
ddist1 <- datadist(BP_women)
options(datadist='ddist1')
k1 <- with(BP_women, quantile(Body.Mass.Index, c(.05, 0.25, .75, .95)))
k1　#　女性で使っていたknot の位置，とその値．

      5%      25%      75%      95% 
19.79504 24.63853 33.28757 41.35564 

idx_model1 <- lrm(Hyt2g_medi_inclu ~ rcs(Body.Mass.Index, k)+Age+Smoking+Ever_drinker+
                   Current_Betel_Chewing + Fruit + Veg,
                 data=BP_women)
#ここは単純にロジスティクス回帰モデルにより当てはまる．
plot(Predict(idx_model1, Body.Mass.Index))

ddist1$limits["Adjust to","Body.Mass.Index"] <- 23
#BMI=23に対する比較するため，reference pointの設定．
idx_model1 <- update(idx_model1)
#単純グラフの場合．
plot(Predict(idx_model1, Body.Mass.Index,ref.zero=TRUE, fun=exp))

dataplot <- Predict(idx_model1,Body.Mass.Index, ref.zero = TRUE, fun=exp)
#ここからは加工したグラフのプログラム：
ggplot(dataplot,aes(Body.Mass.Index, yhat)) +
  #  geom_line(colour="Black", linetype="dashed", size=1.5)+
  theme(plot.subtitle = element_text(vjust = 1),
        plot.caption = element_text(vjust = 1),
        axis.line = element_line(size = 0.5,
                                 linetype = "solid"), panel.grid.major = element_line(colour = "gray98"),
        axis.title = element_text(size = 15),
        panel.background = element_rect(fill = "gray99",
                                        colour = "white", linetype = "twodash"),
        plot.background = element_rect(fill = "white")) +
  scale_x_continuous(breaks = seq(15,50,by=2.5),expression(paste("Body Mass Index", ", ", kg/m^{2})))+
  scale_y_continuous(limits = c(0,10), breaks = seq(0,10,by=0.5),"Odds Ratios (95%CI)")+
  labs(caption = NULL)+
  annotate("text", x=20, y=9, parse = TRUE,
           label="Women",
           size=5)+
  geom_hline(yintercept =1, linetype="dashed")

---

1. Dose-response analyses using restricted cubic spline functions in public health research↩