高維中介數據：基于交替方向乘子法（ADMM）的高維度單模態中介模型的參數估計（入門+實操）

全文摘要

?用于高維度單模態中介模型的參數估計，采用交替方向乘子法（ADMM）進行計算。該包提供了確切獨立篩選（SIS）功能來提高中介效應的敏感性和特異性，并支持Lasso、彈性網絡、路徑Lasso和網絡約束懲罰等不同正則化方法。

Pathway Lasso

背景

傳統的結構方程建模（SEM）在處理大量中介變量時變得不穩定且計算復雜。Pathway Lasso引入了一個新的懲罰函數，它是一種非凸乘積函數的凸松弛，使得同時估計和選擇路徑效應成為可能。通過使用交替方向乘子法（ADMM）的算法，Pathway Lasso可以以閉合形式求解參數，并且其估計器在大樣本下具有漸近一致性。Pathway Lasso的新方法用于在高維中介變量的情況下估計和選擇路徑效應。

實現方法

Pathway Lasso是一種針對高維中介變量問題的新方法，它通過結構方程建模（SEM）的正則化途徑來處理。在高維設置中，當中介變量的數量接近或大于樣本量時，該方法聚焦于估計和選擇路徑效應。為了改善估計的穩定性，Pathway Lasso避免將高維中介變量直接降低為線性組合，這通常是通過主成分分析（PCA）或其他矩陣分解技術實現的，但這些方法限制了對每個中介路徑的解釋性。相反，Pathway Lasso引入了一個新的凸懲罰項，即Pathway Lasso懲罰，直接對路徑效應進行正則化。這種方法解決了傳統Lasso和其他凸正則化方法無法處理的乘積參數問題，因為路徑效應通常表示為兩個參數的乘積，這是一個非凸函數。通過Pathway Lasso懲罰，可以同時實現路徑選擇和路徑效應估計，允許模型直接處理相關中介變量，提供更直接和簡單的中介路徑解釋，尤其適用于分析多個大腦區域作為中介變量的情況。

Pathway Lasso的優勢

在路徑選擇和估計準確性方面相較于其他方法具有以下優勢

高路徑選擇準確性：在模擬數據和fMRI數據集上的應用表明，Pathway Lasso 提出的方法比其他方法具有更高的路徑選擇準確性。
低估計偏誤：Pathway Lasso 方法在估計路徑效應時表現出更低的偏差。
解決非凸性問題：Pathway Lasso 引入了一個新的凸懲罰，直接對乘積非凸函數進行正則化，解決了現有方法未處理的問題。
直接和明確的解釋性：與使用線性組合（如主成分分析）的方法相比，Pathway Lasso 允許對每個中介路徑進行更直接和更簡單的解釋。
處理相關中介變量：Pathway Lasso 允許直接建模相關中介變量，適合分析多個大腦區域作為中介的設置。

實現方法

隨機生成單模態高維度中介分析數據

代碼格式

modalityMediationDataGen(n = 100,p = 50,sigmaY = 1,sizeNonZero = c(3, 3, 4),alphaMean = c(6, 4, 2),alphaSd = 0.1,betaMean = c(6, 4, 2),betaSd = 0.1,sigmaM1 = NULL,gamma = 3,generateLaplacianMatrix = FALSE,seed = 20231201
)

?參數說明

n: 高維中介模型中的主體數量。
p: 高維中介變量的數量。
sigmaY: 因變量誤差分布的標準差。
sizeNonZero: 非零中介變量的數量，生成大、中、小中介效應的模擬場景。
alphaMean, alphaSd: 中介變量與自變量之間效應的平均值和標準差向量。
betaMean, betaSd: 中介變量與因變量之間效應的平均值和標準差向量。
sigmaM1: 中介變量間誤差分布的協方差矩陣，默認為對角矩陣。
gamma: 直接效應的真值。
generateLaplacianMatrix: 邏輯值，指定是否生成網絡懲罰的拉普拉斯矩陣。
seed: 隨機種子，默認為NULL以使用當前種子

返回結果解釋

MediData: 高維中介模型的模擬數據。
MediPara: 中介效應和直接效應的真值。
Info: 輸出包括隨機種子、參數設置以及生成中介模型的拉普拉斯矩陣。

示例代碼?

## 生成分析數據
simuData <- modalityMediationDataGen(seed = 20231201)
str(simuData)
# 輸出結果如下
# List of 3
# $ MediData:List of 3
# ..$ X : num [1:100, 1] 0 0 1 0 0 0 1 0 0 1 ...
# ..$ M1: num [1:100, 1:50] 1.023 -0.369 4.812 1.476 0.188 ...
# ..$ Y : num [1:100, 1] -10.27 6.54 175.08 -1.66 17.55 ...
# $ MediPara:List of 3
# ..$ alpha: num [1, 1:50] 5.99 5.99 6 4.11 4.17 ...
# ..$ beta : num [1:50, 1] 6.11 5.96 6.01 4.05 3.88 ...
# ..$ gamma: num [1, 1] 3
# $ Info    :List of 4
# ..$ parameters     :List of 7
# .. ..$ sigmaY     : num 1
# .. ..$ sizeNonZero: num [1:3] 3 3 4
# .. ..$ alphaMean  : num [1:3] 6 4 2
# .. ..$ alphaSd    : num [1:3] 0.1 0.1 0.1
# .. ..$ betaMean   : num [1:3] 6 4 2
# .. ..$ betaSd     : num [1:3] 0.1 0.1 0.1
# .. ..$ sigmaM1    : num [1:50, 1:50] 1 0 0 0 0 0 0 0 0 0 ...
# ..$ trueValue      :List of 1
# .. ..$ gamma: num [1, 1] 3
# ..$ laplacianMatrix: NULL
# ..$ seed           : num 20231201simuData <- modalityMediationDataGen(seed = 20231201, generateLaplacianMatrix = TRUE)
str(simuData)
simuData <- modalityMediationDataGen(n = 50, p = 1000, seed = 20231201)
str(simuData)

交叉驗證:cvSingleModalityAdmm

通過設置`numFolds`參數進行交叉驗證，可以評估不同懲罰參數下的模型性能，幫助選擇最佳模型

`交叉驗證的結果，用于評估不同參數組合下Pathway Lasso懲罰方法的效果。輸出結果是一個表格，其中包含以下列：

1. **rho**：這是ADMM算法中的ρ參數的候選值，它影響算法的收斂速度和解的質量。
2. **lambda1a**：Pathway Lasso懲罰中的λ1a參數的候選值，L1 范數懲罰中介變量和自變量之間的影響。
3. **lambda1b**：Pathway Lasso懲罰中的λ1b參數的候選值，中介變量和因變量之間影響的 L1 范數懲罰。
4. **lambda1g**：Pathway Lasso懲罰中的λ1g參數的候選值，直接效應的 L1 范數懲罰。默認值為 10 以解決高估問題。
5. **kappa**：Pathway Lasso懲罰的L1范數參數，控制路徑正則化的具體形式。控制了路徑結構的稀疏性，當?kappa?較小時，懲罰的作用更加平滑，有利于保留更多的特征；當?kappa?較大時，懲罰更加集中，有利于稀疏性，即更多特征被剔除。
6. **nu**：Pathway Lasso懲罰的L2范數參數，同樣影響路徑正則化。nu: 控制了路徑結構中特征之間的相關性，當?nu?較小時，路徑結構更加獨立，有利于減少特征之間的相關性；當?nu?較大時，更多的特征將共享相同的路徑，有助于保留相關性較強的特征。
7. **measure**：評估指標，默認均方根誤差（RMSE），用于衡量預測結果與真實結果之間的差異。低的RMSE值通常意味著更好的模型性能，因為這表示預測誤差更小。通過比較這些結果，可以選取最優的參數組合來構建最終模型。

8. lambda2a、lambda2b: 是 Pathway Lasso 方法中額外引入的懲罰項的參數。它們可以控制特征之間的相關性，幫助更好地保留特征間的相關性信息。

lambda2a：L2 范數懲罰中介變量和自變量之間的影響
lambda2b：中介變量和因變量之間影響的 L2 范數懲罰

# 2種不同的懲罰方法## 1.使用交叉驗證進行 ElasticNet 懲罰參數調優
# 執行交叉驗證
cvElasticNetResults <- cvSingleModalityAdmm(X = simuData$MediData$X,  # 獨立變量的數據矩陣（暴露/治療/組）Y = simuData$MediData$Y,  # 因變量的數據向量（結果響應）M1 = simuData$MediData$M1,  # 單模態中介變量numFolds = 5,  # 交叉驗證的折數typeMeasure = "rmse",  # 評估指標類型，默認為均方根誤差rho = c(0.9, 1, 1.1),  # rho 參數的候選值序列lambda1a = c(0.1, 0.5, 1),  # lambda1a 參數的候選值序列lambda1b = c(0.1, 0.3),  # lambda1b 參數的候選值序列lambda1g = c(1, 2),  # lambda1g 參數的候選值序列lambda2a = c(0.5, 1),  # lambda2a 參數的候選值序列lambda2b = c(0.5, 1),  # lambda2b 參數的候選值序列penalty = "ElasticNet"  # 使用 ElasticNet 懲罰
)# 輸出結果： 
> cvElasticNetResultsrho lambda1a lambda1b lambda1g lambda2a lambda2b  measure[1,] 0.9      0.1      0.1        1      0.5      0.5 18.23108[2,] 1.0      0.1      0.1        1      0.5      0.5 18.32964[3,] 1.1      0.1      0.1        1      0.5      0.5 18.17303[4,] 0.9      0.5      0.1        1      0.5      0.5 17.77722[5,] 1.0      0.5      0.1        1      0.5      0.5 17.78040[6,] 1.1      0.5      0.1        1      0.5      0.5 17.77446[7,] 0.9      1.0      0.1        1      0.5      0.5 17.80479
[到達getOption("max.print") -- 略過很多行]]
attr(,"class")
[1] "cvSingleModalityAdmm"--------------------------------------------------------------------------
# 2. 使用交叉驗證進行 Pathway Lasso 懲罰參數調優（lambda2a, lambda2b 未調整）
# 執行交叉驗證
cvPathwayLassoResults <- cvSingleModalityAdmm(X = simuData$MediData$X,  # 獨立變量的數據矩陣（暴露/治療/組）Y = simuData$MediData$Y,  # 因變量的數據向量（結果響應）M1 = simuData$MediData$M1,  # 單模態中介變量numFolds = 5,  # 交叉驗證的折數typeMeasure = "rmse",  # 評估指標類型，默認為均方根誤差rho = c(0.9, 1, 1.1),  # rho 參數的候選值序列lambda1a = c(0.1, 0.5, 1),  # lambda1a 參數的候選值序列lambda1b = c(0.1, 0.3),  # lambda1b 參數的候選值序列lambda1g = c(1, 2),  # lambda1g 參數的候選值序列lambda2a = 1,  # 給定 lambda2a 參數值lambda2b = 1,  # 給定 lambda2b 參數值penalty = "PathwayLasso",  # 使用 Pathway Lasso 懲罰penaltyParameterList = list(kappa = c(0.5, 1), nu = c(1, 2))  # 懲罰參數列表，包括 kappa 和 nu
)# 輸出結果：
cvPathwayLassoResultsrho lambda1a lambda1b lambda1g kappa nu  measure[1,] 0.9      0.1      0.1        1   0.5  1 19.46943[2,] 1.0      0.1      0.1        1   0.5  1 19.37725[3,] 1.1      0.1      0.1        1   0.5  1 19.40920[4,] 0.9      0.5      0.1        1   0.5  1 19.49747
[到達getOption("max.print") -- 略過很多行]]
attr(,"class")
[1] "cvSingleModalityAdmm"

將權矩陣轉換為拉普拉斯矩陣的輔助函數:weightToLaplacian()?

# 將權矩陣轉換為拉普拉斯矩陣的輔助函數:weightToLaplacian() 
set.seed(20231201) # 設置隨機數種子
p <- 5 # 設置節點數
W <- matrix(0, nrow = p, ncol = p) # 初始化權矩陣
W[lower.tri(W)] <- runif(p*(p-1)/2, 0, 1) # 生成隨機權的下三角矩陣
W[upper.tri(W)] <- t(W)[upper.tri(W)] # 使權矩陣對稱
diag(W) <- 1 # 對角線元素設為1
W
# 輸出結果如下
# [,1]      [,2]       [,3]      [,4]       [,5]
# [1,] 1.0000000 0.1623753 0.48119340 0.4406640 0.36219565
# [2,] 0.1623753 1.0000000 0.41138920 0.1344408 0.64471664
# [3,] 0.4811934 0.4113892 1.00000000 0.5306324 0.08042435
# [4,] 0.4406640 0.1344408 0.53063239 1.0000000 0.85450197
# [5,] 0.3621956 0.6447166 0.08042435 0.8545020 1.00000000(L <- weightToLaplacian(W)) # 將權矩陣轉換為拉普拉斯矩陣
# 輸出結果如下
# [,1]        [,2]        [,3]        [,4]        [,5]
# [1,]  0.59124083 -0.06767837 -0.19443191 -0.16374871 -0.13501050
# [2,] -0.06767837  0.57499652 -0.16949748 -0.05094059 -0.24505056
# [3,] -0.19443191 -0.16949748  0.60058145 -0.19491464 -0.02963414
# [4,] -0.16374871 -0.05094059 -0.19491464  0.66218945 -0.28956112
# [5,] -0.13501050 -0.24505056 -0.02963414 -0.28956112  0.66007653

擬合高維單模態中介模型

根據cvSingleModalityAdmm的結果挑選最佳參數，擬合🔤高維單模態中介模型🔤

penalty方法

penalty方法有3種+ 各自對應的懲罰參數列表【penaltyParameterList】

默認為彈性網絡 ElasticNet
lambda1a, lambda1b, lambda1g, lambda2a, lambda2b
路徑套索（PathywayLasso）
kappa 路徑 Lasso 的 L1 范數懲罰。
nu 路徑 Lasso 的 L2 范數懲罰
網絡約束懲罰（Network）
需要應用于網絡懲罰的拉普拉斯矩陣

確定獨立性篩選?(SIS)

SIS：指定是否執行確定獨立性篩選 (sure independence screening, SIS)

SISThreshold,中介者目標降維的閾值。默認值為“2”，這會將維度減少到 2*n/log(n)。n代表樣本量

輸出結果

gamma：🔤估計直接影響🔤
alpha：🔤估計中介變量和自變量之間的影響。🔤
beta：🔤估計中介變量和因變量之間的影響🔤

綜合應用

1.?ElasticNet 懲罰

## 生成經驗數據
simuData <- modalityMediationDataGen(seed = 20231201, generateLaplacianMatrix = TRUE)## ElasticNet 懲罰的參數估計
modelElasticNet <- singleModalityAdmm( X = simuData$MediData$X, Y = simuData$MediData$Y, M1 = simuData$MediData$M1, rho = 1, lambda1a = 1, lambda1b = 0.1, lambda1g = 2, lambda2a = 1, lambda2b = 1, penalty = "ElasticNet" )# 擬合并預測
fitted(modelElasticNet) 
predict(modelElasticNet, matrix(c(0, 1), ncol=1))# SIS獨立性篩選
simuData <- modalityMediationDataGen(n = 50, p = 1000, seed = 20231201)
modelElasticNetSIS <- singleModalityAdmm( X = simuData$MediData$X, Y = simuData$MediData$Y, M1 = simuData$MediData$M1, rho = 1, lambda1a = 1, lambda1b = 0.1, lambda1g = 2, lambda2a = 1, lambda2b = 1, penalty = "ElasticNet", SIS = TRUE ) fitted(modelElasticNetSIS) 
predict(modelElasticNetSIS, matrix(c(0, 1), ncol=1))

2.?使用拉普拉斯矩陣進行網絡懲罰的參數估計


# 1.使用模擬數據中的拉普拉斯矩陣
simuData <- modalityMediationDataGen(seed = 20231201, generateLaplacianMatrix = TRUE)modelNetwork <- singleModalityAdmm( X = simuData$MediData$X, Y = simuData$MediData$Y, M1 = simuData$MediData$M1, rho = 1, lambda1a = 1, lambda1b = 0.1, lambda1g = 2, lambda2a = 1, lambda2b = 1, penalty = "Network", penaltyParameterList = list(laplacianMatrix = simuData$Info$laplacianMatrix) )# 2. 自定義的拉普拉斯矩陣set.seed(20231201) 
p <- ncol(simuData$MediData$M1) 
W <- matrix(0, nrow = p, ncol = p) 
W[lower.tri(W)] <- runif(p*(p-1)/2, 0, 1) 
W[upper.tri(W)] <- t(W)[upper.tri(W)] 
diag(W) <- 1 
L <- weightToLaplacian(W) modelNetwork <- singleModalityAdmm( X = simuData$MediData$X, Y = simuData$MediData$Y, M1 = simuData$MediData$M1, rho = 1, lambda1a = 1, lambda1b = 0.1, lambda1g = 2, lambda2a = 1, lambda2b = 1, penalty = "Network", penaltyParameterList = list(laplacianMatrix = L) )

3.?Pathway Lasso 懲罰的參數估計

simuData <- modalityMediationDataGen(seed = 20231201, generateLaplacianMatrix = TRUE)modelPathwayLasso <- singleModalityAdmm( X = simuData$MediData$X, Y = simuData$MediData$Y, M1 = simuData$MediData$M1, rho = 1, lambda1a = 1, lambda1b = 0.1, lambda1g = 2, lambda2a = 1, lambda2b = 1, penalty = "PathwayLasso", penaltyParameterList = list(kappa = 1, nu = 2) )