R语言 PLS Path Modeling 模型诊断及验证
我们最后完成路径模型的简单构建，我们假设潜变量之间的关系是，越高的攻击质量和越好的防守质量，能带来更多胜利。模型初步得出的结论与假设的结论不一致。
完成模型的简单构建仅仅是第一步，接下来要做的更加重要：检查获得的结果并评估模型的质量。本篇文章讨论如何解释结果及诊断模型。在此之前，我们先初步了解plspm()函数，这对我们后续调整并诊断模型非常重要。
plspm(Data, path_matrix, blocks, modes = NULL,scaling = NULL, scheme = "centroid", scaled = TRUE,tol = 1e-06, maxiter = 100, plscomp = NULL,boot.val = FALSE, br = NULL, dataset = TRUE)
'plspm'包中的plspm() 函数有13个参数。可以根据目的将其分为三个主要部分，第一部分与定义PLS路径模型的参数有关，第二部分与PLS算法有关，第三部分是其他选项。可使用help()函数查看函数的详细信息。
help("plspm")定义PLS路径模型
Data ##数据，包含度量数据的数组或矩阵
path_matrix ##定义内部模型的矩阵
blocks ##定义外部模型变量区块的列表
scaling ##定义非度量数据变量的度量范围的列表
modes ##定义每个区块的测量模式的向量
PLS算法
scheme ##内部模型的加权方案
scaled ##是否应将数据标准化
tol ##用于检查PLS迭代阶段的收敛算法，公差阈值
maxiter ##最大迭代次数，默认为100
plscomp ##当处理非度量数据时，表示PLS组件的数量其他选项
boot.val ##指示是否必须执行引导验证
br ##引导程序重采样次数
dataset ##是否应检索数据矩阵
大多数情况下，只需调整其中几个参数。最重要的部分是前三个参数：Data, path_matrix, blocks。
了解plspm()函数的基础使用后，运行代码得到路径模型的初步结果并尝试解读。
# load package plspmlibrary(plspm)# load data spainfootdata(spainfoot)# rows of the path matrixAttack = c(0, 0, 0)Defense = c(0, 0, 0)Success = c(1, 1, 0)# inner model matrixfoot_path = rbind(Attack, Defense, Success)# add column namescolnames(foot_path) = rownames(foot_path)# blocks of indicators (outer model)foot_blocks = list(1:4, 5:8, 9:12)# vector of modes (reflective)foot_modes = rep("A", 3)# run plspm analysisfoot_pls = plspm(spainfoot, foot_path, foot_blocks, modes = foot_modes)
模型评估
PLS路径模型有内部模型和外部模型组成，模型的评估需要对内部模型和外部模型进行分析和解释。诊断遵循两个步骤：
（1）对外部模型进行评估
（2）对内部模型进行评估
查看结果：
# what's in foot_pls?foot_pls## Partial Least Squares Path Modeling (PLS-PM)## ---------------------------------------------## NAME DESCRIPTION## 1 $outer_model outer model## 2 $inner_model inner model## 3 $path_coefs path coefficients matrix## 4 $scores latent variable scores## 5 $crossloadings cross-loadings## 6 $inner_summary summary inner model## 7 $effects total effects## 8 $unidim unidimensionality## 9 $gof goodness-of-fit## 10 $boot bootstrap results## 11 $data data matrix## ---------------------------------------------## You can also use the function 'summary'# summarized resultssummary(foot_pls)
外部模型评估
对外部模型进行评估首先应该区分每个潜变量的度量指标是形成（or formative indicators）还是反映（reflective indicators ）,前者的评估较为简单，如果存在多重共线性，应对指标进行选择。后者的评估比较复杂，需要遵循三个方面：指标单一性（unidimensionality ）；指标由其潜在变量正确解释；评估给定结构与其他结构的差异程度。接下来将分别从这三个方面对模型进行诊断。
unidimensionality
指标单一性，plspm()构建模型得出结果后，利用 $unidim 查看
# unidimensionalityunidim# Mode MVs C.alpha DG.rho eig.1st eig.2nd# Attack A 4 0.8906 0.92456 3.017 0.7923# Defense A 4 0.0000 0.02602 2.393 1.1753# Success A 4 0.9165 0.94233 3.217 0.5370#第一列表示测量类型，第二列表示每个区块中MVs的数量#第三列 C.alpha 第四列 DG.rho 以及后两列的eig.1st eig.2nd 用于衡量单一性
一般认为，cronbach's alpha>0.7时可以接受。在我们的例子中，Attack和Success的C.alpha值大于0.7，Defense则不满足这个条件，说明用于衡量Defense的manifest variables可能有潜在的错误以至于不满足指标的单一性原则。另一种衡量区块指标单一性的度量标准是Dillon-Goldstein’s rho，同样的，当DG.rho值大于0.7时认为区块指标是满足指标单一性原则，Defense区块的DG.rho值为0.02602，也反应了Defense区块的指标有问题。第三个度量指标涉及每个区块指标相关矩阵的特征值，分别用到了第一个特征值（eig.1st）和第二个特征值（eig.2nd），如果区块满足指标单一性原则，其第一个特征值应该大于1，第二个特征值应该小于1，Defense区块不满足这个条件。
# cronbach's alphaunidim[, 3, drop = FALSE]# C.alpha# Attack 0.8906# Defense 0.0000# Success 0.9165# dillon-goldstein rhounidim[, 4, drop = FALSE]# DG.rho# Attack 0.92456# Defense 0.02602# Success 0.94233# eigenvaluesunidim[, 5:6]# eig.1st eig.2nd# Attack 3.017 0.7923# Defense 2.393 1.1753# Success 3.217 0.5370
从上面的分析可以了解到Defense区块存在问题，我们可以通过plot()函数并设置参数 what="loadings" 来观察外部模型。
# plotting loadingsplot(foot_pls, what = "loadings")
从图中可以看到。Defense区块指标中，CSA和CSH箭头为红色，表明他们与Defense区块的关系是负相关，同样的我们可以检查外部模型以确认。
## outer model resultsfoot_pls$outer_model## name block weight loading communality redundancy## 1 GSH Attack 0.3366 0.9379 0.8797 0.0000## 2 GSA Attack 0.2819 0.8621 0.7432 0.0000## 3 SSH Attack 0.2892 0.8408 0.7070 0.0000## 4 SSA Attack 0.2396 0.8263 0.6828 0.0000## 5 GCH Defense -0.1087 0.4837 0.2339 0.0000## 6 GCA Defense -0.3915 0.8759 0.7672 0.0000## 7 CSH Defense 0.3274 -0.7464 0.5571 0.0000## 8 CSA Defense 0.4035 -0.8926 0.7968 0.0000## 9 WMH Success 0.2309 0.7755 0.6014 0.5145## 10 WMA Success 0.3030 0.8864 0.7856 0.6722## 11 LWR Success 0.2821 0.9686 0.9382 0.8027## 12 LRWL Success 0.2958 0.9437 0.8906 0.7620## 使用subset()函数可以帮助我们只观察外部模型中defense部分的结果## Defense outer model resultssubset(foot_pls$outer_model, block == "Defense")## name block weight loading communality redundancy## 5 GCH Defense -0.1087 0.4837 0.2339 0.0000## 6 GCA Defense -0.3915 0.8759 0.7672 0.0000## 7 CSH Defense 0.3274 -0.7464 0.5571 0.0000## 8 CSA Defense 0.4035 -0.8926 0.7968 0.0000
可以观察到权重（weights）和loading，通过plot()函数可视化。
# plotting weightsplot(foot_pls, what = "weights")
可以看到Defense权重图中，一半指标的权重为正，一半指标的权重为负，这就是其cronbach's alpha和Goldstein’s rho值较小的原因。此外，这也是导致defense与success之间的路径系数为负的原因。这是矛盾的，意味着越好的防御质量，获胜就越少。
尽管GCH（主场失球总数）和GCA（客场失球总数）与防守有关，但是其衡量的是防守的不足。越高的GCH和GCA，意味着防守质量越差。我们需要调整指标。一种方法是在原有的数据中添加新的两列，NGCH和NGCA，计算为负的GCH和负的GCA，调整后的指标与的defense区块就有了正相关关系，并使用调整后的指标来重新构建路径分析模型。
# add two more columns NGCH and NGCAspainfoot$NGCH = -1 * spainfoot$GCHspainfoot$NGCA = -1 * spainfoot$GCA# check column namesnames(spainfoot)## [1] "GSH" "GSA" "SSH" "SSA" "GCH" "GCA" "CSH" "CSA" "WMH"## [10] "WMA" "LWR" "LRWL" "YC" "RC" "NGCH" "NGCA# new list of blocks (with column positions of variables)new_blocks_pos = list(1:4, c(15,16,7,8), 9:12)# new list of blocks (with names of variables)new_blocks_str = list(c("GSH", "GSA", "SSH", "SSA"),c("NGCH", "NGCA", "CSH", "CSA"),c("WMH", "WMA", "LWR", "LRWL"))# re-apply plspmfoot_pls = plspm(spainfoot, foot_path, new_blocks_str, modes = foot_modes)# plot loadingsplot(foot_pls, "loadings")
调整指标后的外部模型：
重新检查模型指标单一性：
## unidimensionalityfoot_pls$unidim## Mode MVs C.alpha DG.rho eig.1st eig.2nd## Attack A 4 0.8906 0.9246 3.017 0.7923## Defense A 4 0.7718 0.8549 2.393 1.1753## Success A 4 0.9165 0.9423 3.217 0.5370
Loadings and Communalities
接下来，检查外部模型的Loadings 和 Communalities，Loadings表示潜在变量与其指标的相关性，Communalities代表的是相关的平方，目的是检查区块中潜在变量是否能够很好的解释指标。
## loadings and communalitiesfoot_pls$outer_model## name block weight loading communality redundancy## 1 GSH Attack 0.3366 0.9380 0.8798 0.0000## 2 GSA Attack 0.2819 0.8621 0.7432 0.0000## 3 SSH Attack 0.2893 0.8408 0.7070 0.0000## 4 SSA Attack 0.2396 0.8263 0.6828 0.0000## 5 NGCH Defense 0.1088 0.4837 0.2339 0.0000## 6 NGCA Defense 0.3914 0.8759 0.7672 0.0000## 7 CSH Defense 0.3274 0.7464 0.5571 0.0000## 8 CSA Defense 0.4035 0.8926 0.7967 0.0000## 9 WMH Success 0.2309 0.7755 0.6014 0.5145## 10 WMA Success 0.3030 0.8864 0.7856 0.6722## 11 LWR Success 0.2821 0.9686 0.9382 0.8027## 12 LRWL Success 0.2958 0.9437 0.8906 0.7620
良好的指标相关性（Loading）应大于0.7。在示例中，Defense区块中变量NGCH的相关性为0.4837，小于0.7。可以考虑删去指标。在实际应用中，如果遇到相关性不佳的指标，应根据情况选择保留或删除对应的指标。
Cross-loadings
除了检查指标与其对应的潜在变量之间的相关性，还必须检查所谓的Cross-loadings，即指标与其他潜在变量之间的相关性。例如，Attack区块中GSH指标，其相关性为0.9380，这个值应当大于后两个值。如果一个指标与其他潜在变量之间的相关性大于其与对应潜在变量的相关性，需要考虑这个指标的适当性。
# cross-loadingsfoot_pls$crossloadings## name block Attack Defense Success## 1 GSH Attack 0.9380 0.5159 0.8977## 2 GSA Attack 0.8621 0.3391 0.7519## 3 SSH Attack 0.8408 0.4139 0.7714## 4 SSA Attack 0.8263 0.3361 0.6390## 5 NGCH Defense 0.1305 0.4837 0.1598## 6 NGCA Defense 0.4622 0.8759 0.5751## 7 CSH Defense 0.3188 0.7464 0.4810## 8 CSA Defense 0.4215 0.8926 0.5928## 9 WMH Success 0.7086 0.4226 0.7755## 10 WMA Success 0.7731 0.7115 0.8864## 11 LWR Success 0.8444 0.5380 0.9686## 12 LRWL Success 0.8601 0.5892 0.9437
也可以通过ggplot2包进行可视化。
# load ggplot2 and reshapelibrary(ggplot2)library(reshape)# reshape crossloadings data.frame for ggplotxloads = melt(foot_pls$crossloadings, id.vars = c("name", "block"),variable_name = "LV")# bar-charts of crossloadings by blockggplot(data = xloads,aes(x = name, y = value, fill = block)) +# add horizontal reference linesgeom_hline(yintercept = 0, color = "gray75") +geom_hline(yintercept = 0.5, color = "gray70", linetype = 2) +# indicate the use of car-chartsgeom_bar(stat = 'identity', position = 'dodge') +# panel display (i.e. faceting)facet_wrap(block ~ LV) +# tweaking some grahical elementstheme(axis.text.x = element_text(angle = 90),line = element_blank(),plot.title = element_text(size = 12)) +# add titleggtitle("Crossloadings")
可视化：
内部模型结构评估
在评估外部模型后，进行内部模型结构评估。检查结构方程中每个回归的结果，查看结果使用 $inner_model，示例中只有一个回归元素，即success
# inner modelinner_model# $Success# Estimate Std. Error t value Pr(>|t|)# Intercept -2.356e-16 0.09217 -2.556e-15 1.000e+00# Attack 7.573e-01 0.10440 7.254e+00 1.349e-06# Defense 2.836e-01 0.10440 2.717e+00 1.466e-02
此外，还要通过其他三个指标来评估结构模型的质量：R2，冗余指标和GOF。
R2
用R2评估模型质量：示例的R2为0.855，根据标准表明模型良好。
低：R < 0.30
中：0.30 < R < 0.60
高：R > 0.60
## inner model summaryinner_summary# Type R2 Block_Communality Mean_Redundancy AVE# Attack Exogenous 0.0000 0.7532 0.0000 0.7532# Defense Exogenous 0.0000 0.5887 0.0000 0.5887# Success Endogenous 0.8556 0.8040 0.6878 0.8040select R2inner_summary[, "R2", drop = FALSE]# R2# Attack 0.0000# Defense 0.0000# Success 0.855
冗余指标
高冗余指标意味着高预测能力，示例中success的Mean_Redundancy为0.6878，,表示attack和defense能够解释success指标变异性的68%。
GOF
GOF指标是为拟合优度度量，同时考虑了测量和结构模型的模型质量。良好模型的GOF值应大于0.7。示例中GoF值为0.78可以解释为该模型的预测能力为78％。
# gof indexgof# [1] 0.7823
验证结果
由于PLS-PM不基于任何分布假设，因此参数估计的显著性水平（基于正态分布）不适用于PLS-PM。取而代之的是，使用重采样来获取有关参数估计值可变性的信息。
自举（bootstrapping）是一种非参数方法，用于估计PLS参数估计的精度。简单描述过程可以理解为：创建M个样本，以获得PLS模型中每个参数的M个估计。每个样本都是通过从原始数据集中进行替换采样而获得的，样本大小等于原始数据集中的案例数。
通过plsm()函数中的 boot.val 来控制，函数默认重采样次数为100，假如我们将重采样次数设置为200，代码如下：
# running bootstrap validationfoot_val = plspm(spainfoot, foot_path, new_blocks_str, modes = foot_modes,boot.val = TRUE, br = 200)# bootstrap resultsfoot_val$boot
得到的结果中，分别是参数的原始值，平均值，标准误差以及置信区间为95%的较低百分比和较高百分比。
以上就是PLS-PM模型构建的结果的诊断和验证。
最后将调整后的路径模型绘制出来：
plot(foot_pls)参考文献：
http://www.gastonsanchez.com/PLS Path Modeling withR.pdf