マルチレベル分析がなぜ必要かの説明をRで簡単に実装して行う
臨床をやっている時には意識していなかったが、施設差を解析に入れないと解析結果が明らかにmisleadingとなる事があるので、Rで実装しておく
個人的には、「その治療法は施設によってやり方が全然違うんちゃうん？」もしくは「その病気の予後は地域によって全然違うんちゃうん？」と思ったら、マルチレベルはやるべきだと思う。
マルチレベル分析とは？
マルチレベル分析とは施設差を考慮する際によく使われる解析方法である。過去にある人の研究発表を聞いた後に「施設差を考慮しないといけないのでは？」と質問した際に「いや、これだけの症例数あるから大丈夫でしょ」みたいな事を言われた事があるのだが、この記事を読んで(Rで実験して)症例数を増やすだけでは全く大丈夫ではない事が伝わればと思う。
まずはデータベース作成のため、以下をコピペする。
library(mvtnorm) N <- 1000; se1 <- 150; se2 <- -90 sigma = matrix(c(se1, se2, se2, se1), ncol = 2) rand = rmvnorm(n = N, mean = c(130, 360), sigma) x0 <- rep("A",N) x1 = rand[, 1] x2 = rand[, 2] X1 <- cbind(data.frame(x0),x1,x2) rand_2 = rmvnorm(n = N, mean = c(140, 380), sigma) x0 <- rep("B",N) x1 = rand_2[, 1] x2 = rand_2[, 2] X2 <- cbind(data.frame(x0),x1,x2) rand_3 = rmvnorm(n = N, mean = c(120, 340), sigma) x0 <- rep("C",N) x1 = rand_3[, 1] x2 = rand_3[, 2] X3 <- cbind(data.frame(x0),x1,x2) X_total <- rbind(X1,X2,X3) colnames(X_total) <- c("state","sBP","time")
作成したデータベースの説明と解析目的
解析目的：「食事指導にかけた時間」と「血圧」の関係を調べたいとする1
データベースには
str(X_total) #'data.frame': 3000 obs. of 3 variables: # $ state: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ... # $ sBP : num 149 141 111 122 123 ... # $ time : num 341 359 372 378 358 ... head(X_total,3) # state sBP time #1 A 148.8035 341.0151 #2 A 140.5801 358.7447 #3 A 110.9342 371.7364
- “state”は3つの地域の名称が
- “sBP”には収縮期血圧が
- “time”には年あたりの減塩にかけた指導時間が
入っているとする
stateのイメージとしては、Aは普通の地域、Bは医療リテラシーが低い(不健康な)地域、Cはリテラシーが高い地域とする
実際解析してみる
施設を無視して(多施設をまぜこぜにして)血圧と指導時間を単変量解析し、散布図に直線を引く
res_total <- lm(sBP~time,data=X_total) summary(res_total) ## ## Call: ## lm(formula = sBP ~ time, data = X_total) ## ## Residuals: ## Min 1Q Median 3Q Max ## -58.428 -9.577 0.041 10.122 48.306 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 90.94149 4.64591 19.575 <2e-16 *** ## time 0.10935 0.01288 8.487 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 14.44 on 2998 degrees of freedom ## Multiple R-squared: 0.02346, Adjusted R-squared: 0.02314 ## F-statistic: 72.03 on 1 and 2998 DF, p-value: < 2.2e-16 plot(X_total[,"time"],X_total[,"sBP"]) abline(res_total,col=2)
timeの係数が有意差をもって正になったので、 食事指導にかけた時間が長ければ長いほど、血圧が上がってしまうという事が検証されてしまった。(plotでも傾きは正となっている)
これは施設の違いを考慮せずに解析した事が影響している。
施設を考慮すると？
施設毎の性質(の分の分散)は別のパラメータに引き取ってもらうので、きちんと減塩指導時間の効果が解析できるというイメージ。
まず最初に素データを施設毎に色分けして散布図をみると、
plot(X_total[,"time"],X_total[,"sBP"],col=as.numeric(as.factor(X_total[,"state"])))
散布図を見るだけでもなんとなく違うんじゃないかと思えてくるが、ここに施設ごとに回帰式を作って線を引いてみる
d_a <- subset(X_total,X_total$state %in% "A") res_a <- lm(sBP~time,data=d_a) d_b <- subset(X_total,X_total$state %in% "B") res_b <- lm(sBP~time,data=d_b) d_c <- subset(X_total,X_total$state %in% "C") res_c <- lm(sBP~time,data=d_c) plot(X_total[,"time"],X_total[,"sBP"],col=as.numeric(as.factor(X_total[,"state"]))) c="magenta" abline(res_a,col=c) abline(res_b,col=c) abline(res_c,col=c)
施設毎にみると、指導時間が長い方が血圧は下がるという結果になる。
初めのスクリプトで仮にN=1000000としても変わらないので対象患者数を増やせば良いという話ではないという事はわかるだろう。
解析結果として出力してみる
まずは施設を考慮せずに解析した結果は、
summary(res_total) ## ## Call: ## lm(formula = sBP ~ time, data = X_total) ## ## Residuals: ## Min 1Q Median 3Q Max ## -58.428 -9.577 0.041 10.122 48.306 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 90.94149 4.64591 19.575 <2e-16 *** ## time 0.10935 0.01288 8.487 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 14.44 on 2998 degrees of freedom ## Multiple R-squared: 0.02346, Adjusted R-squared: 0.02314 ## F-statistic: 72.03 on 1 and 2998 DF, p-value: < 2.2e-16
となり、有意差をもって指導時間が増えると血圧が高いという結果になる。
次に施設差を加えてマルチレベル分析をすると、
summary(lme4::lmer(sBP~time+(1|state), data=X_total)) ## Linear mixed model fit by REML ['lmerMod'] ## Formula: sBP ~ time + (1 | state) ## Data: X_total ## ## REML criterion at convergence: 22103.1 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -3.8264 -0.6823 0.0100 0.6831 3.6717 ## ## Random effects: ## Groups Name Variance Std.Dev. ## state (Intercept) 490.23 22.141 ## Residual 92.02 9.593 ## Number of obs: 3000, groups: state, 3 ## ## Fixed effects: ## Estimate Std. Error t value ## (Intercept) 346.15689 13.78917 25.10 ## time -0.59956 0.01435 -41.77 ## ## Correlation of Fixed Effects: ## (Intr) ## time -0.375
となり、初めに設定した共分散の値に近づいたのがわかる。2
マルチレベルで考慮するのは、地域変数・病院変数・外科の手術なら術者毎のカテゴリ・内科なら外来毎のカテゴリなどが挙げられるが、臨床やってる頃は知らんかった。これはRCTをやる時もランダム性を保てない要因に関しては同様に必要だと思う。
