投稿者: zanet_infozanet

応答変数（従属変数）が二値変数、例えばある疾患・そうでないで、いくつかの説明変数からそれを推定するモデルを作成したい場合、ロジスティック回帰分析を用いることができます。説明変数は連続変数でも名義変数でも使えます。

応答変数が3つ以上ある場合はMulti-nomial logistic regression analysis多項ロジスティック回帰分析を用いることができます。例えば、診断したい病態が3つあり、各症例で3種類の検査のデータがあるような場合、このデータを解析して、各検査の値に対する係数を算出し、別の症例や検証用の症例のデータに対して、その症例に対する疾患確率を計算することができます。診断したい疾患が3つであれば、それぞれに対する疾患確率の値を3つ得ることができます。実臨床では、ある疾患・そうでないの2つを鑑別することよりも、3つ以上の複数の疾患のどれかを診断するような状況の方がはるかに多いと思います。前者は感度・特異度のパラダイムで処理できますが、後者の場合、感度・特異度の考え方では対処が難しくなります。

想定される疾患・病態が3つ以上あり、検査や診断法が2つ以上あるような場合、Multi-nomial logistic regressionに基づくモデルは役に立つはずです。

Rのパッケージでnnetを用いると、その解析ができます。解析に必要なデータは図1のようなものになります。

このようなデータをExcelなどで用意しておきます。なお、ここに示すのは架空のデータです。

nnetはBrian Ripley, William VenablesによるRのパッケージです。 Link

以下Rでスクリプトを実行していきます。変数名は図1のものを用いていますが、実際に図1のデータを解析しているわけではありません。スクリプトと手順を示すのが目的です。それぞれ自分のデータで必要に応じて変数名など変更して実行してください。

#Install packages needed.

packneed=c(“nnet”,”tcltk2″);current=installed.packages();addpack=setdiff(packneed,rownames(current));url=”https://cran.ism.ac.jp/”;if(length(addpack)>0){install.packages(addpack,repos=url)};if(length(addpack)==0){print(“Already installed.”)}

今回取り上げたパッケージ以外に、NeuralNetTools、EffectStars2というパッケージも有用です。今回Multi-nomial logistic regression analysisをnnetパッケージのmultinom()関数で行っていますが、nnetはそれだけでなく、Neural networkを用いる解析もできます。

#Read in packages.

library(nnet)
library(tcltk2)

#Read in data.

train.dat=read.delim(“clipboard”,sep=”\t”,header=TRUE)

#Excelで必要なデータ範囲をコピーしてからRに戻りこのスクリプトを実行します。クリップボードを介して、データが変数train.datに格納されます。Macの場合は、train.dat=read.delim(pipe(“pbpaste”),sep=”\t”,header=TRUE) です。

#Set a reference as a factor variable.

train.dat$Diagnosis = as.factor(train.dat$Diagnosis)
train.dat$Diagnosis = relevel(train.dat$Diagnosis, ref = “D1”)

ベースとなるDiagnosisを設定します。その際に、数値ではなくfactorとして取り扱うので、as.factor()関数を用いています。ref引数で設定するDiagnosisは通常は対照となるような診断・病態です。

#Multi-nomial logistic regression analysis.

multinom_model = multinom(Diagnosis ~ Test1+Test2+Test3, data = train.dat)

診断 ~ 説明変数1 + 説明変数2 + 説明変数3の様にformulaを記述します。

#Show the results: Intercepts, coefficients for explanatory variables.

summary(multinom_model)

#P values for explanatory variables.

z=summary(multinom_model)$coefficients/summary(multinom_model)$standard.errors
P_value=(1-pnorm(abs(z),0,1))*2
P_value

summary()関数ではP値は計算されないので、このようなスクリプトで計算しています。

#Odds ratios for explanatory variables (exponentials of the coefficients).

exp(coef(multinom_model))

#Prediction the diagnosis with explanatory variables for each case with the training data.

(predDiagnosis=predict(multinom_model, newdata=train.dat, type=”probs”))

ここでは、トレーニングに用いたデータに対して、Predictionを行っています。適合度が高いので、正診率は高くなります。

#Prepare the result data to be saved.

results=cbind(train.dat,predDiagnosis)

#Save the result data as a text file (tab-separated).

filnam=tclvalue(tkgetSaveFile(initialfile=”.txt”,filetypes=”{{Text Files} {.txt}} {{All files} *}”));if(filnam!=””){write.table(results,filnam,sep=”\t”,row.names=FALSE,col.names=TRUE)}

ファイル名を設定して任意のフォルダーにテキストファイルとして、元のデータとPredictionの結果を合わせてタブ区切りのテーブルの形式で保存します。

#Calcuating disease probabilities for individual explanatory variable data.
Set data for the case:
an example

Test1=3
Test2=14
Test3=5
case=cbind(Test1,Test2,Test3)
row.names(case)[1]=”case”

任意の検査結果の場合の、診断の推定を行います。

#Output Diagnosis class.

(CasePredicted=predict(multinom_model, newdata = case, type=”class”))

3つの診断の内、疾患確率が一番高いものが出力されます。

#Output probabilities.

(CasePredicted=predict(multinom_model, newdata = case, type=”probs”))

3つの診断のそれぞれの疾患確率が出力されます。

診断のPredictionは上記のpredict()関数を用いてできますが、自分でスクリプトを書いて計算したい場合は以下の様にします。

Multi-nomial logistic regressionの結果、Reference diagnosis以外のDiagnosisに対する切片、係数の値が得られます。これらの値を用いて、各症例の疾患確率を計算します。切片、係数の値は、coef(multinom_model)でデータフレームとして得られます。そこから、必要な切片、係数の値を読み出して、その症例のそれぞれの検査結果の値にそれぞれの対応する係数を掛け算して、総和を求め、さらに切片の値を加算します。得られた値のExponentialを求めます。それを各Diagnosisについて計算し、それらの総和を求め、さらに1を加算します。その値で、各DiagnosisのExponentialの値を割り算すると、各Diagnosisに対する疾患確率が得られます。加算する1はReference diagnosisの切片と係数を0に設定するためで、0のExponentialが1になるためです。

検査が3つ、すなわち、説明変数が3つの場合で、その症例の検査結果がx1, x2, x3でああれば、
診断iについて、yi = ai + b1i*x1 + b2i*x2 + b3i*x3
P(診断i) = exp(yi)/(1 + Σexp(yi))

通常のLogistic regressionの場合は、y = a + b1*x1 + b2*x2 + b3*x3で、
P(診断)=1/(1 + exp(-y)) ですから計算法が違います。

文献：
Ciaburro G, Venkateswaran B: Neural Networks with R: Smart models using CNN, RNN, deep learning, and artificial intelligence principles. 2017, Packt Publishing, Birmingham, UK.

Xu J: Modern Applied Regression. 2023, CRC press, Taylor & Francis Group.

Multi-nomial logistic regression with R. R-bloggers.