LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@GORIN

【R コードの説明】コヨーテのマイクロバイオーム解析 ~ SEM（構造化方程式モデリング）~

Posted at

カナダにおいて、16S rRNA遺伝子アンプリコンシークエンシングを用いて、糞便微生物叢が都市部のコヨーテにおける人為的な食餌（ファストフードの廃棄物など）と不健康の関係、また同様に農村部のコヨーテにおけるタンパク質が豊富な食餌（小型草食動物や偶蹄類など）と健康の関係をどのように媒介するかを調べた研究を見つけました。

<今回> 本研究の最後では、SEM（構造方程式モデリング）を用いて、食餌と微生物分類の関係をモデル化している。今回は、本研究で実行されたSEM部分の各コードの説明を付記していきます。

※本研究で構築されたモデルによって、以下の関係が支持されました。
「都市部のコヨーテがより多くの人為的な食餌(ファストフードなどの廃棄物)を摂取し、それが多様性の増加と腸球菌および連鎖球菌の存在量の増加を通じて、健康に間接的な悪影響を及ぼす。」(2020, Sugden et al.)

Rスクリプトおよび論文のリンク

SEM(Structural Equation Modeling)

#### STRUCTURAL EQUATION MODELS ####
# (1) Prepare data for SEM ####
feces.SEM.clr <- cbind(subset(comp.strict.clr[["Genus"]], is.na(Observed)=="FALSE"), feces.strict.sum.clr[["Family"]])

# Convert factors to dummy variables for lavaan.
levels(feces.SEM.clr$Em.PCR.Overall)[levels(feces.SEM.clr$Em.PCR.Overall)=="Y"] <- 1
levels(feces.SEM.clr$Em.PCR.Overall)[levels(feces.SEM.clr$Em.PCR.Overall)=="N"] <- 0
feces.SEM.clr$Em.PCR.Overall <- as.numeric(as.character(feces.SEM.clr$Em.PCR.Overall))

levels(feces.SEM.clr$Location)[levels(feces.SEM.clr$Location)=="urban"] <- 1
levels(feces.SEM.clr$Location)[levels(feces.SEM.clr$Location)=="rural"] <- 0
feces.SEM.clr$Location <- as.numeric(as.character(feces.SEM.clr$Location))

levels(feces.SEM.clr$Sex)[levels(feces.SEM.clr$Sex)=="M"] <- 1
levels(feces.SEM.clr$Sex)[levels(feces.SEM.clr$Sex)=="F"] <- 0
feces.SEM.clr$Sex <- as.numeric(as.character(feces.SEM.clr$Sex))

# Scale stomach content variables to meet a normal distribution.
feces.SEM.clr$Vol.Prey <- log(feces.SEM.clr$Vol.Prey+0.01)
feces.SEM.clr$Vol.Anthro <- log(feces.SEM.clr$Vol.Anthro+0.01)

# Reduce the magnitude of some variables so that variances are equal.
feces.SEM.clr$Vol.Anthro <- feces.SEM.clr$Vol.Anthro/10
feces.SEM.clr$Vol.Prey <- feces.SEM.clr$Vol.Prey/10
feces.SEM.clr$Observed_Extrap <- feces.SEM.clr$Observed_Extrap/100
feces.SEM.clr$PD_Extrap <- feces.SEM.clr$PD_Extrap/10

# (2) SEM FOR ANTHROPOGENIC FOOD / POOR HEALTH CONNECTIONS ####
# For simplicity, only the initial model and the final model are shown.
# From the null model, additional paths were added as recommended by modification indices.
# Non-significant (p>0.1) paths were removed.
# All models were evaluated based on AIC, NNFI, SRMR, RMSEA, and CFI, as shown.
# Model selection stopped when the addition or removal of a path caused model AIC to increase,
# or other fit parameters (NNFI, SRMR, RMSEA, CFI) to get worse, 
# even if the path being removed was not significant.
# Criteria for NNFI, SRMR, RMSEA, CFI, Chi-sq: https://www.cscu.cornell.edu/news/Handouts/SEM_fit.pdf 

# Create storage areas for model texts and model fits.
div.gen.models <- list()
div.gen.fits <- list()

# Construct the initial model.
div.gen.models[[1]] <- '
# regressions
index1.PC1 ~ Location + d13C + Vol.Anthro + Vol.Anthro + Em.PCR.Overall + Enterococcus + Streptococcus + PD_Extrap
Spleen.to.Body.Ratio ~ Location + d13C + Vol.Anthro + Em.PCR.Overall + PD_Extrap + Enterococcus + Streptococcus
d13C ~ Location
Vol.Anthro ~ Location
PD_Extrap ~ Vol.Anthro + d13C + Location
Enterococcus ~ Location + d13C + Vol.Anthro
Streptococcus ~ Location + d13C + Vol.Anthro
Em.PCR.Overall ~ Location + d13C + Vol.Anthro + PD_Extrap
Enterococcus ~~ Streptococcus
PD_Extrap ~~ Streptococcus
PD_Extrap ~~ Enterococcus'

# Calculate and evaluate the SEM fit.
div.gen.fits[[1]] <- sem(div.gen.models[[1]], data=feces.SEM.clr)
summary(div.gen.fits[[1]], standardized=TRUE)

# Determine if any paths should be added.
mi <- modindices(div.gen.fits[[1]])
print(mi[mi$mi > 3.0,])

# Run the next model after adding or removing paths, based on the previous model results.
div.gen.models[[2]] <- '
# regressions
index1.PC1 ~ Streptococcus + PD_Extrap
Spleen.to.Body.Ratio ~ Location + Em.PCR.Overall + Enterococcus
d13C ~ Location
PD_Extrap ~ d13C + Location + Vol.Anthro
Enterococcus ~ Vol.Anthro
Streptococcus ~ Vol.Anthro
Em.PCR.Overall ~ PD_Extrap
Enterococcus ~~ Streptococcus
PD_Extrap ~~ Streptococcus
PD_Extrap ~~ Enterococcus'

div.gen.fits[[2]] <- sem(div.gen.models[[2]], data=feces.SEM.clr)
summary(div.gen.fits[[2]], standardized=TRUE)
mi <- modindices(div.gen.fits[[2]])
print(mi[mi$mi > 3.0,])

# Compare the two models.
anova(div.gen.fits[[2]], div.gen.fits[[1]])

# Evaluate model AIC.
aictab.lavaan(div.gen.fits,
              paste0("Model", seq(length(div.gen.fits))))


            
# (3) SEM FOR PROTEIN / GOOD HEALTH CONNECTIONS ####
# Create storage areas for model texts and model fits.
protein.gen.models <- list()
protein.gen.fits <- list()

# Propose the initial model.
protein.gen.models[[1]] <- '
# regressions
index1.PC1 ~ Location + d15N + Vol.Prey + Vol.Prey + Em.PCR.Overall + Anaerobiospirillum + Fusobacterium + Sutterella
Spleen.to.Body.Ratio ~ Location + d15N + Em.PCR.Overall + Anaerobiospirillum + Fusobacterium + Sutterella
d15N ~ Location
Vol.Prey ~ Location
Anaerobiospirillum ~ Location + d15N + Vol.Prey
Fusobacterium ~ Location + d15N + Vol.Prey
Sutterella ~ Location + d15N + Vol.Prey
Em.PCR.Overall ~ Location + d15N + Vol.Prey
Anaerobiospirillum ~~ Fusobacterium
Anaerobiospirillum ~~ Sutterella
Sutterella ~~ Fusobacterium'

# Calculate and evaluate the SEM fit.
protein.gen.fits[[1]] <- sem(protein.gen.models[[1]], data=feces.SEM.clr)
summary(protein.gen.fits[[1]], standardized=TRUE)

# Determine if any paths should be added.
mi <- modindices(protein.gen.fits[[1]])
print(mi[mi$mi > 3.0,])

# Run the next model.
# (Accelerated here to show the final model).
protein.gen.models[[2]] <- '
# regressions
index1.PC1 ~ Location + d15N
Spleen.to.Body.Ratio ~ Location + Em.PCR.Overall + Anaerobiospirillum + Sutterella
d15N ~ Location
Vol.Prey ~ Spleen.to.Body.Ratio
Anaerobiospirillum ~ d15N
Fusobacterium ~ Location + d15N
Sutterella ~ d15N
Em.PCR.Overall ~ Location + Vol.Prey
Anaerobiospirillum ~~ Fusobacterium
Anaerobiospirillum ~~ Sutterella
Sutterella ~~ Fusobacterium'

protein.gen.fits[[2]] <- sem(protein.gen.models[[2]], data=feces.SEM.clr)
summary(protein.gen.fits[[2]], standardized=TRUE)
mi <- modindices(protein.gen.fits[[2]])
print(mi[mi$mi > 3.0,])

# Compare to the previous model.
anova(protein.gen.fits[[6]], protein.gen.fits[[5]])

# Evaluate model AICs.
aictab.lavaan(protein.gen.fits,
              paste0("Model", seq(length(protein.gen.fits))))

# Evaluate additional fit statistics.
protein.gen.fits.stats <- data.frame(matrix(nrow=42))
for(i in 1:length(protein.gen.fits)){
  vector <- fitmeasures(protein.gen.fits[[i]])
  protein.gen.fits.stats[,i] <- vector
  rownames(protein.gen.fits.stats) <- names(vector)
}
protein.gen.fits.stats <- subset(protein.gen.fits.stats, rownames(protein.gen.fits.stats) %in% c("npar", "chisq", "df", "pvalue",
                                                                                                 "cfi", "nnfi", "aic", "rmsea", "srmr", "gfi"))
colnames(protein.gen.fits.stats) <- paste0("Model", seq(length(protein.gen.fits)))
for(i in 1:ncol(protein.gen.fits.stats)){protein.gen.fits.stats[,i] <- as.numeric(protein.gen.fits.stats[,i])}
protein.gen.fits.stats <- as.data.frame(t(protein.gen.fits.stats))

# Determine the R-sq value for each variable.
lavInspect(protein.gen.fits[[6]], "rsquare")

# Export to CytoScape for figure design.
t <- semPlot::semPaths(protein.gen.fits[[6]], "std", title=FALSE, layout="spring", 
                       residuals=FALSE, intercepts=FALSE, thresholds=FALSE, nCharNodes=0, sizeMan=10)
t <- igraph::as.igraph(t)
igraph::write.graph(t, file = "~/protein.gen.SEM.gml", format = "gml")

# (4) SEM FOR SPLEEN MASS CONNECTIONS ####
spleen.models <- list()
spleen.fits <- list()

spleen.models[[1]] <- '
# regressions
index1.PC1 ~ Location + d13C + Vol.Anthro + Em.PCR.Overall + Erysipelotrichaceae + Coriobacteriaceae + Lachnospiraceae
Spleen.to.Body.Ratio ~ Location + d13C + Vol.Anthro + Em.PCR.Overall + Erysipelotrichaceae + Coriobacteriaceae + Lachnospiraceae
d13C ~ Location
Vol.Anthro ~ Location
Erysipelotrichaceae ~ Location + d13C + Vol.Anthro
Coriobacteriaceae ~ Location + d13C + Vol.Anthro
Lachnospiraceae ~ Location + d13C + Vol.Anthro
Em.PCR.Overall ~ Location + d13C + Vol.Anthro
Erysipelotrichaceae ~~ Coriobacteriaceae
Coriobacteriaceae ~~ Lachnospiraceae
Erysipelotrichaceae ~~ Lachnospiraceae'

spleen.fits[[1]] <- sem(spleen.models[[1]], data=feces.SEM.clr)
summary(spleen.fits[[1]], standardized=TRUE)
mi <- modindices(spleen.fits[[1]])
print(mi[mi$mi > 3.0,])

spleen.models[[9]] <- '
# regressions
index1.PC1 ~ Location + Lachnospiraceae
Spleen.to.Body.Ratio ~ Location + Em.PCR.Overall + Erysipelotrichaceae
d13C ~ Location
Erysipelotrichaceae ~ Location + Vol.Anthro
Coriobacteriaceae ~ Location + Vol.Anthro
Lachnospiraceae ~ Location
Em.PCR.Overall ~ Location + d13C
Erysipelotrichaceae ~~ Coriobacteriaceae
Coriobacteriaceae ~~ Lachnospiraceae
Erysipelotrichaceae ~~ Lachnospiraceae'

spleen.fits[[9]] <- sem(spleen.models[[9]], data=feces.SEM.clr)
summary(spleen.fits[[9]], standardized=TRUE)
mi <- modindices(spleen.fits[[9]])
print(mi[mi$mi > 3.0,])

anova(spleen.fits[[9]], spleen.fits[[1]])

aictab.lavaan(spleen.fits,
              paste0("Model", seq(length(spleen.fits))))

spleen.fits.stats <- data.frame(matrix(nrow=42))
for(i in 1:length(spleen.fits)){
  vector <- fitmeasures(spleen.fits[[i]])
  spleen.fits.stats[,i] <- vector
  rownames(spleen.fits.stats) <- names(vector)
}
spleen.fits.stats <- subset(spleen.fits.stats, rownames(spleen.fits.stats) %in% c("npar", "chisq", "df", "pvalue",
                                                                                  "cfi", "nnfi", "aic", "rmsea", "srmr", "gfi"))
colnames(spleen.fits.stats) <- paste0("Model", seq(length(spleen.fits)))
for(i in 1:ncol(spleen.fits.stats)){spleen.fits.stats[,i] <- as.numeric(spleen.fits.stats[,i])}
spleen.fits.stats <- as.data.frame(t(spleen.fits.stats))

lavInspect(spleen.fits[[9]], "rsquare")

t <- semPlot::semPaths(spleen.fits[[9]], "std", title=FALSE, layout="spring", 
                       residuals=FALSE, intercepts=FALSE, thresholds=FALSE, nCharNodes=0, sizeMan=10)
t <- igraph::as.igraph(t)
igraph::write.graph(t, file = "~/spleen.SEM.gml", format = "gml")

コードの詳細

  • (1) Prepare data for SEM
    • この部分では、SEMを行うためのデータの準備をしています。
    • feces.SEM.clrというデータフレームを作成しています。このデータフレームは、comp.strict.clrというデータフレームのGenusという列と、feces.strict.sum.clrというデータフレームのFamilyという列を結合したものです。また、Observedという列に欠損値がある行は除外しています。

  • Covert factors to dummy variables for lavaan
    • feces.SEM.clrのEm.PCR.Overall, Location, Sexという列は因子型の変数ですが、lavaanパッケージでSEMを行うためには、ダミー変数に変換する必要があります。そのため、levels関数とas.numeric関数を使って、それぞれの列の値を0と1に置き換えています。例えば、Em.PCR.Overallという列では、Yという値を1に、Nという値を0に変換します。

  • Scala stomach content variables to meet a normal distribution.
    • feces.SEM.clrのVol.PreyとVol.Anthroの列は、胃の内容物の量を表す変数ですが、正規分布に従っていない可能性があります。そのため、log関数を使って、対数変換を行っています。また、0.01を足していますが、これは0の値を持つ行がある場合に、log(0)が定義できないのを防ぐためです。

  • Reduce the magnitude of some variable so that variance are equal.
    • feces.SEM.clrのVol.Anthro, Vol.Prey, Observed_Extrap, PD_Extrapという列は、変数の大きさが他の変数と比べて大きすぎると分散が不均一になる可能性があります。そのため、それぞれの列の値を10や100で割って、変数の大きさを小さくしています。

  • (2) SEM FOR ANTHROPOGENIC FOOD / POOR HEALTH CONNECTIONS
    • この部分では、人為的な食物と健康状態の関係を調べるためのSEMを行っています。

  • Create storage areas for model texts and model fits.
    • div.gen.modelsとdiv.gen.fitsというリストを作成しています。これらのリストは、モデルのテキストとモデルの適合度を格納するためのものです。

  • Construct the initial model.
    • div.gen.models[[1]]という文字列に、初期モデルの式を記述しています。このモデルでは、index1.PC1とSpleen.to.Body.Ratioという2つの潜在変数を仮定しています。index1.PC1は、腸内細菌の多様性を表す変数で、Genusという列の主成分分析の第1主成分です。Spleen.to.Body.Ratioは、脾臓の大きさと体重の比率で、免疫系の状態を表す変数です。これらの潜在変数は、Location, d13C, Vol.Anthro, Em.PCR.Overall, Enterococcus, Streptococcus, PD_Extrapという観測変数に影響を与えると仮定しています。Locationは都市か田舎かを表す変数で、d13Cは炭素同位体比で、食物の起源を表す変数です。Vol.Anthroは人為的な食物の量で、Vol.Preyは天然の獲物の量です。Em.PCR.Overallは、エンテロバクター科の細菌の有無を表す変数で、EnterococcusとStreptococcusは、それぞれエンテロコッカス属とストレプトコッカス属の細菌の量を表す変数です。PD_Extrapは、腸内細菌の多様性の推定値です。このモデルでは、観測変数間の回帰式と、潜在変数間の共分散を記述しています。

  • Calculate and evaluate the SEM fit.
    • div.gen.fits[[1]]というオブジェクトに、sem関数を使って、初期モデルの適合度を計算しています。sem関数の引数には、モデルの式とデータフレームを指定しています。
    • summary関数を使って、初期モデルの適合度の要約を表示します。standardized=TRUEというオプションを使って、標準化されたパラメータ推定値も表示します。

  • Determine if any paths should be added.
    • modindices関数を使って、初期モデルの修正指数を計算します。修正指数とは、モデルにパスを追加したり削除したりした場合に、どの程度モデルの適合度が改善されるかを示す指標です。修正指数が大きいほど、モデルの改善の余地があることを意味します。print関数を使って、修正指数が3.0より大きいパスを表示しています。
    • div.gen.models[[2]]という文字列に、修正指数に基づいて、初期モデルからパスを追加したり削除したりした最終モデルの式を記述しています。このモデルでは、index1.PC1 と Spleen.to.Body.Ratio という2つの潜在変数を仮定していますが、観測変数に与える影響が変わっています。index1.PC1は、StreptococcusとPD_Extrapに影響を与えると仮定しています。Spleen.to.Body.Ratioは、Location, Em.PCR.Overall, Enterococcusに影響を与えると仮定しています。d13Cは、LocationとPD_Extrapに影響を与えると仮定しています。Vol.Anthroは、PD_Extrap, Enterococcus, Streptococcusに影響を与えると仮定しています。Em.PCR.Overallは、PD_Extrapに影響を与えると仮定しています。このモデルでは、観測変数間の回帰式と、潜在変数間の共分散を記述しています。
    • div.gen.fits[[2]]というオブジェクトに、sem関数を使って、最終モデルの適合度を計算しています。sem関数の引数には、モデルの式とデータフレームを指定しています。
    • summary関数を使って、最終モデルの適合度の要約を表示しています。standardized=TRUEというオプションを使って、標準化されたパラメータ推定値も表示しています。
    • modindices関数を使って、最終モデルの修正指数を計算しています。修正指数が3.0より大きい場合を表示していますが、このモデルでは、修正指数が3.0より大きいパスはありません。これは、モデルの改善の余地がないことを意味します。

  • Compare the two models
    • anova関数を使って、初期モデルと最終モデルを比較しています。anova関数は、モデル間の適合度の差を検定するために、カイ二乗統計量の差分と自由度の差分を用います。anova関数の結果から、最終モデルは初期モデルよりも有意に適合度が高いことを確かめます。

  • Evaluate model AIC.

    • aictab.lavaan関数は、モデル間の適合度をAICという指標で比較します。AICは、モデルの複雑さと当てはまりの良さをバランスさせるための情報量基準です。AICが小さいほど、モデルの適合度が高いと判断されます。aictab.lavaan関数の結果から、最終モデルは初期モデルよりもAICが小さいこと確かめます。
    • 以上から、最終モデルは初期モデルよりも優れたモデルであると言えます。最終モデルでは、潜在変数のindex1.PC1とSpleen.to.Body.Ratioが観測変数に与える影響の大きさや方向が明らかになります。また、観測変数間の相関や潜在変数間の相関も推定されます。※
      これらのパラメータ推定値は、summary関数の出力やsemPlot関数の図で確認できます。semPlot関数は、構造方程式モデルをグラフィカルに表示するための関数です。semPlot関数の引数には、モデルの適合度オブジェクトと、ラベルの種類、レイアウト、スタイル、ノードの文字数、ノードのサイズ、エッジのラベルのサイズなどを指定できます。

  • Evaluate additional fit statistics.

    • この部分では、protein.gen.fitsというリストに格納されたモデルの適合度を評価するためのコードを書いています。
    • protein.gen.fits.statsというデータフレームを作成しています。このデータフレームは、42行の空の行列から作られています。
    • for文を使って、protein.gen.fitsの各要素に対して、fitmeasures関数を適用しています。fitmeasures関数は、モデルの適合度の指標を計算するための関数です。fitmeasures関数の結果をvectorという変数に格納しています。vectorの名前をprotein.gen.fits.statsの行名に、vectorの値をprotein.gen.fits.statsの列に代入しています。protein.gen.fits.statsの列名は、Model1, Model2, …というように番号をつけています。
    • protein.gen.fits.statsの行をsubset関数を使って、必要な指標だけに絞り込んでいます。必要な指標は、npar（パラメータの数）、chisq（カイ二乗統計量）、df（自由度）、pvalue（p値）、cfi（比較適合指数）、nnfi（非規範適合指数）、aic（赤池情報量基準）、rmsea（平方根平均二乗誤差近似）、srmr（標準化残差平均二乗）、gfi（適合指数）です。
    • protein.gen.fits.statsの各要素をas.numeric関数を使って、数値に変換しています。
    • protein.gen.fits.statsをt関数を使って、転置しています。これは、行と列を入れ替えるためです。

  • Determine the R-sq value for each variable.
    • この部分では、protein.gen.fits[[6]]というオブジェクトに格納されたモデルの、各観測変数の決定係数（R-sq）を計算するためのコードを書いています。
    • lavInspect関数を使って、protein.gen.fits[[6]]というオブジェクトの"rsquare"という属性を取得しています。lavInspect関数は、モデルの適合度オブジェクトの中身を調べるための関数です。"rsquare"という属性は、各観測変数の決定係数を表しています。決定係数とは、観測変数の分散のうち、モデルで説明できる割合を示す指標です。決定係数が高いほど、モデルの当てはまりが良いと判断されます。

  • Export to CytoScape for figure design.
    • この部分では、protein.gen.fits[[6]]というオブジェクトに格納されたモデルの図を作成して、CytoScapeというソフトウェアにエクスポートするためのコードを書いています。
    • semPlot::semPaths関数を使って、protein.gen.fits[[6]]というオブジェクトのモデルをグラフィカルに表示しています。semPlot::semPaths関数は、構造方程式モデルを図で表すための関数です。semPlot::semPaths関数の引数には、モデルの適合度オブジェクトと、標準化されたパラメータ推定値を使うかどうか（“std”）、タイトルを表示するかどうか（title=FALSE）、レイアウトの種類（layout=“spring”）、残差や切断点などを表示するかどうか（residuals=FALSE, intercepts=FALSE, thresholds=FALSE）、ノードの文字数（nCharNodes=0）、ノードのサイズ（sizeMan=10）などを指定しています。
    • semPlot::semPaths関数の結果をtという変数に格納しています。tという変数は、構造方程式モデルの図を表すオブジェクトです。
    • igraph::as.igraph関数を使って、tという変数をigraphというパッケージで扱えるオブジェクトに変換しています。igraphというパッケージは、ネットワーク分析を行うためのパッケージです。igraph::as.igraph関数は、様々な形式のネットワークデータをigraphオブジェクトに変換するための関数です。
    • igraph::write.graph関数を使って、igraphオブジェクトをファイルに書き出しています。igraph::write.graph関数は、ネットワークデータを様々な形式で保存するための関数です。igraph::write.graph関数の引数には、ネットワークデータと、ファイル名と、ファイル形式を指定しています。ここでは、"~/protein.gen.SEM.gml"というファイル名と、"gml"というファイル形式を指定しています。gmlというファイル形式は、CytoScapeというソフトウェアで読み込める形式です。CytoScapeというソフトウェアは、ネットワークデータの可視化や分析を行うためのソフトウェアです。

  • (3) SEM FOR PROTEIN / GOOD HEALTH CONNECTIONS
    ここでは、タンパク質と健康の関係を調べるためSEMを実行します。

  • Create storage areas for model texts and model fits.

    • モデルのテキストとフィットを保存するためのリストを作成します。

  • Propose the initial model.

    • 初期モデルを作成します。このモデルでは、index1.PC1とSpleen.to.Body.Ratioという二つの従属変数に対して、Location, d15N, Vol.Prey, Em.PCR.Overall, Anaerobiospirillum, Fusobacterium, Sutterellaという七つの独立変数と、Anaerobiospirillum, Fusobacterium, Sutterellaという三つの潜在変数の関係を回帰式で表現しています。また、潜在変数同士の共分散も指定しています。

  • Calculate and evaluate the SEM fit

    • sem関数を用いて、初期モデルのフィットを計算し、summary関数で標準化されたパラメータ推定値やモデル適合度指標を表示します。

  • Determine if any paths should be added.

    • modindices関数を用いて、モデル改善のために追加すべきパスを探します。モディフィケーションインデックス（MI）が3.0以上のパスを表示します。

  • Run the next model.

    • モデルを実行します。このモデルでは、初期モデルから不要なパスを削除し、Vol.PreyをSpleen.to.Body.Ratioの関数として表現しています。これは、MIの値に基づいて行われています。
    • sem関数とsummary関数で、次のモデルのフィットとパラメータ推定値を得ます。また、modindices関数で、さらに改善できるパスを確認します。

  • Compare to the previous model.

    • anova関数で、初期モデルと次のモデルを比較します。モデル間のカイ二乗差検定を行い、有意水準を判断します。

  • Evaluate model AICs.

    • aictab.lavaan関数で、モデルのAICを評価します。AICは、モデルの複雑さと適合度のバランスをとる指標で、小さいほど良いモデルとされます。

  • Evaluate additional fit statistics.

    • fitmeasures関数で、モデルの他の適合度指標を評価します。nparはモデルのパラメータ数、chisqはモデルのカイ二乗値、dfは自由度、pvalueはモデルのp値、cfiは比較適合指数、nnfiは非常適合指数、rmseaは平方根平均二乗誤差近似、srmrは標準化残差平均二乗根、gfiは適合指数です。これらの指標は、モデルのデータへの当てはまりの良さを表します。

  • Determine the R-sq value for each variables.

    • lavInspect関数で、各変数の決定係数（R-sq）を求めます。決定係数は、変数の分散がモデルでどれだけ説明されているかを示す指標です。

  • Export to CytoScape for figure design.

    • semPlot::semPaths関数で、モデルのパス図を作成します。この図は、変数間の因果関係やパラメータ推定値を視覚的に表現します。
    • igraph::as.igraph関数とigraph::write.graph関数で、パス図をgml形式のファイルとして出力します。このファイルは、CytoScapeというソフトウェアで図のデザインを行うために使用できます。

  • (4) SEM FOR SPLEEN MASS CONNECTIONS
    • この部分では、spleen.modelsとspleen.fitsというリストに格納されたモデルのコードを書いています。これらのモデルは、脾臓の質量とその他の変数の関係を調べるためのものです。
    • spleen.models[[1]]という文字列に、初期モデルの式を記述しています。このモデルでは、index1.PC1とSpleen.to.Body.Ratioという2つの潜在変数を仮定しています。index1.PC1は、腸内細菌の多様性を表す変数で、Genusという列の主成分分析の第1主成分です。Spleen.to.Body.Ratioは、脾臓の大きさと体重の比率で、免疫系の状態を表す変数です。これらの潜在変数は、Location, d13C, Vol.Anthro, Em.PCR.Overall, Erysipelotrichaceae, Coriobacteriaceae, Lachnospiraceaeという観測変数に影響を与えると仮定しています。Locationは都市か田舎かを表す変数で、d13Cは炭素同位体比で、食物の起源を表す変数です。Vol.Anthroは人為的な食物の量で、Vol.Preyは天然の獲物の量です。Em.PCR.Overallは、エンテロバクター科の細菌の有無を表す変数で、Erysipelotrichaceae, Coriobacteriaceae, Lachnospiraceaeは、それぞれエリシペロトリキス科、コリオバクテリウム科、ラクノスピラセウス科の細菌の量を表す変数です。このモデルでは、観測変数間の回帰式と、潜在変数間の共分散を記述しています。
    • spleen.fits[[1]]というオブジェクトに、sem関数を使って、初期モデルの適合度を計算しています。sem関数の引数には、モデルの式とデータフレームを指定しています。
    • summary関数を使って、初期モデルの適合度の要約を表示しています。standardized=TRUEというオプションを使って、標準化されたパラメータ推定値も表示しています。
    • modindices関数を使って、初期モデルの修正指数を計算しています。修正指数とは、モデルにパスを追加したり削除したりした場合に、どの程度モデルの適合度が改善されるかを示す指標です。修正指数が大きいほど、モデルの改善の余地があることを意味します。print関数を使って、修正指数が3.0より大きいパスを表示しています。
    • spleen.models[[9]]という文字列に、修正指数に基づいて、初期モデルからパスを追加したり削除したりした最終モデルの式を記述しています。このモデルでは、index1.PC1とSpleen.to.Body.Ratioという2つの潜在変数を仮定していますが、観測変数に与える影響が変わっています。index1.PC1は、LocationとLachnospiraceaeに影響を与えると仮定しています。Spleen.to.Body.Ratioは、Location, Em.PCR.Overall, Erysipelotrichaceaeに影響を与えると仮定しています。d13Cは、Locationに影響を与えると仮定しています。Vol.Anthroは、ErysipelotrichaceaeとCoriobacteriaceaeに影響を与えると仮定しています。Em.PCR.Overallは、Locationとd13Cに影響を与えると仮定しています。このモデルでは、観測変数間の回帰式と、潜在変数間の共分散を記述しています。
    • spleen.fits[[9]]というオブジェクトに、sem関数を使って、最終モデルの適合度を計算しています。sem関数の引数には、モデルの式とデータフレームを指定しています。
    • summary関数を使って、最終モデルの適合度の要約を表示しています。standardized=TRUEというオプションを使って、標準化されたパラメータ推定値も表示しています。
    • modindices関数を使って、最終モデルの修正指数を計算しています。修正指数が3.0より大きいパスを表示していますが、このモデルでは、修正指数が3.0より大きいパスはありません。これは、モデルの改善の余地がないことを意味します。
    • anova関数を使って、初期モデルと最終モデルを比較しています。anova関数は、モデル間の適合度の差を検定するために、カイ二乗統計量の差分と自由度の差分を用います。anova関数の結果から、最終モデルは初期モデルよりも有意に適合度が高いことを確かめます。
    • aictab.lavaan関数は、モデル間の適合度をAICという指標で比較します。AICは、モデルの複雑さと当てはまりの良さをバランスさせるための基準です。AICが小さいほど、モデルの適合度が高いと判断されます。aictab.lavaan関数から、最終モデルが初期モデルよりもAICが小さいことを確かめます。
    • 以上から、最終モデルが初期モデルよりも優れたモデルであることを言います。最終モデルでは、潜在変数のindex1.PC1とSpleen.to.Body.Ratioが観測変数に与える影響の大きさや方向が明らかになります。また、観測変数間の相関や潜在変数間の相関も推定されます。これらのパラメータ推定値は、summary関数の出力やsemPlot関数の図で確認できます。
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?