Multicollinearity & Heteroskedasticity

Large N & Leeuwenhoek (70700173)

Yue Hu

Overview

1. Collinearity
   • Violation
   • Diagnosis
   • Adjustment
2. Heteroscedasticity
   • Violation
   • Diagnosis
   • Adjustment

• Hetero-: From the Greek word “heteros,” meaning “different” or “other.”
• Skedasticity: Derived from the Greek word “skedastikos,” meaning “able to scatter(spread).”
• It refers to the dispersion or variability of the error term in a regression model. So, heteroskedasticity literally means “different scattering” or “variable dispersion,” which implies that the variance of the error term is not constant.

1 Multicollinearity

1.1 Before going into details

Why should we use multiple regressions, and why controls?

One observes a set of covariates that is, after statistical adjustment, sufficient to make treatment status as-if random.—Keele et al. 2020

1.2 Statistical Adjustment

1. Linearity in the parameter; ← ⍻ Specification, Moderation & GLM
2. Nonstochastic X (“given X,” a.k.a., “X is fixed”); ← ✓ Data collection
3. X has positive noninfinite variance (var(X)); ← ✓ Data collection
4. Correct specification; ← ✓ Specification
5. Identification (N > K; K = 2 for a simple OLS); ← ✓ Data collection & perfect collinearity (Today)
6. Mean zero errors (E(ε_i|X_i) = 0); ← ⍻ Specification (Omited variables, etc.) & Missing
7. No covariance between X_i and ε_i; ← ✗ Endogeneity issue
8. No autocorrelation; ← ✗ Autocorrelation issue
9. No perfect collinearity ← ✗ (Today)
10. Homoskedasticity; ← ✗ (Today)

1.3 Perfect Collinearity

\[X_{2i} = \hat\delta_0 + \hat\delta_1X_{1i} + \hat r_{1i}.\]

When cov(X₁, X₂) = 1, r_1i = 0, then \(\hat\beta_1 = \frac{\sum\hat r_{1i}y_i}{\hat r_{1i}^2}\) cannot be estimated.

E.g., \(X_{2i} = 1 + 2X_{1i}.\)

\[\begin{align} Y_i =& \hat\beta_0 + \hat\beta_1X_{1} + \hat\beta_2X_{2} + \hat \epsilon, \text{(PRF)}\\ =& \hat\beta_0 + \hat\beta_1X_{1} + \hat\beta_2(1 + 2X_1) + \hat \epsilon, \\ =& (\hat\beta_0 + \hat\beta_2) + (\hat\beta_1 + 2\hat\beta_2)X_1 + \hat \epsilon,\\ \Rightarrow Y_i =& \tilde\beta_0 + \tilde\beta_1X_{1} + \hat \epsilon. \end{align}\]

• Consequence
  • \(\hat\epsilon\) is unbiased
  • PRF (β₂) is nonidentifiable

PRF: population regression function

1.4 Multicollinearity

• Coefficient, \(\hat\beta = (\boldsymbol{X'X})^{-1}\boldsymbol{X'Y}\)
  • Small changes in X can lead to dramatic changes in \((\boldsymbol{X'X})^{-1}\), and thus \(\hat\beta\).

Example

\[X = \begin{bmatrix} 1 & 0.99 \\ 0.99 & 1 \end{bmatrix}, \tilde X = \begin{bmatrix} 1 & 0.99 \\ 0.99 & 1.01 \end{bmatrix} \Rightarrow X^{-1} \approx \begin{bmatrix} 50.25 & -49.75 \\ -49.75 & 50.25 \end{bmatrix}, \tilde X^{-1} \approx \begin{bmatrix} 33.44 & 0.3311 \\ 0.3311 & 33.78 \end{bmatrix}\]

• Variance
  • \(var(\hat\beta_1|X) = \sigma^2(\boldsymbol{X'X})^{-1}= \frac{\sigma^2}{\sum(X_{1i} - \bar X_1)^2(1 - \beta_{12}^2)},\) where \(\beta_{12} = cov(X_1, X_2)\).
  • As β₁₂ increases, the variance also increases (Means?)

E.g., \(X = \begin{bmatrix} 1 & 0.99 \\ 0.99 & 1 \end{bmatrix}, \tilde X = \begin{bmatrix} 1 & 0.99 \\ 0.99 & 1.01 \end{bmatrix}\)

To calculate the inverse matrix of matrix X, \(X^{-1} = \frac{1}{\text{det}(X)} \cdot \text{adj}(X)\), where det(X) is the determinant of matrix X, and the adj(X) is the adjugate of matrix X (the transpose of the cofactor matrix of X).

• \(\text{det}(X) = \sum_{i=1}^{n} (-1)^{i+j} \cdot x_{ij} \cdot \text{det}(X_{ij})\) where:
  • \(n\) is the order (number of rows or columns) of matrix X.
  • \(x_{ij}\) represents the element in the \(i\)-th row and \(j\)-th column of matrix X.
  • \(\text{det}(X_{ij})\) is the determinant of the submatrix obtained by removing the \(i\)-th row and \(j\)-th column from X.
  • For a 2×2 matrix, \(\text{det}(X) = x_{11} \cdot x_{22} - x_{12} \cdot x_{21}\)
• \(\text{adj}(X) = (C_{ij})^{T}\), \(C_{ij} = (-1)^{i+j} \cdot \text{det}(X_{ji})\), where:
  • \(C_{ij}\) is the element in the \(i\)-th row and \(j\)-th column of the cofactor matrix \(C\).
  • \(\text{det}(X_{ji})\) is the determinant of the submatrix obtained by removing the \(j\)-th row and \(i\)-th column from X (note the reversal of indices compared to the determinant formula).
• For the X in the example,
  • \(\text{det}(X) = 1 \cdot 1 - 0.99 \cdot 0.99 = 1 - 0.9801 = 0.0199,\)
  • \(\text{adj}(X) = \begin{bmatrix} 1 & -0.99 \\ -0.99 & 1 \end{bmatrix}.\)
  • \(X^{-1} = \frac{1}{0.0199} \cdot \begin{bmatrix} 1 & -0.99 \\ -0.99 & 1\end{bmatrix} \approx \begin{bmatrix} 50.25 & -49.75 \\ -49.75 & 50.25 \end{bmatrix}\)
• For the \(\tilde X\) in the example,
  • \(\text{det}(\tilde X) = 1 \cdot 1.01 - 0.99 \cdot 0.99 = 1.01 - 0.9801 = 0.0299,\)
  • \(\text{adj}(\tilde X) = \begin{bmatrix} 1 & -0.99 \\ -0.99 & 1.01 \end{bmatrix}.\)
  • \(\tilde X^{-1} = \frac{1}{0.0299} \cdot \begin{bmatrix} 1 & -0.99 \\ -0.99 & 1.01 \end{bmatrix} \approx \approx \begin{bmatrix} 33.44 & 0.3311 \\ 0.3311 & 33.78 \end{bmatrix}\)

1.5 Illustration with simulated data

E.g., \(X'_2 = 0.5X_1 + 0.5X_2 + \nu.\)

n <- 10000

df_sim <- tibble(
  x1 = rnorm(n, mean = 3.13), x2 = rnorm(n, sd = 0.4),
  y = 2 + .2 * x1 + .4 * x2 + rnorm(n)
)

df_fake <- df_sim |> 
  mutate(x2 = 0.5 * x1 + 0.5 * x2 + rnorm(n, sd = 0.1))  

cor(df_fake$x1, df_fake$x2)

[1] 0.9118169

1.6 Consequence of multicollinearity: Illustration

ls_result <- map(list(df_sim, df_fake), \(dataset){
  lm(y ~ x1 + x2, data = dataset)
}) |>
  set_names(c("True", "W. multicol"))

small_multiple(ls_result)

1.7 Diagnosis

Variance Inflation Factors (VIF, [1, +∞]): A measure of how much the variance of the estimated coefficient β_x is “inflated” by the correlation among the predictor variables.

\[VIF = \frac{1}{1 - {R}^{2}_{X_i}} = \frac{1}{Tolerance}\]

• \({R}^{2}_{X_i}\): The coefficient of determination for a regression model where X_i is the dependent variable, and all other predictor variables are independent variables.
  • 1: no correlation between X₁ and the rest
  • Rule of thumb: 4, 10
  • E.g., VIF of the true model is 1.0000473, 1.0000473 in the previous slide, and for the one with multicollinearity is 5.9315546, 5.9315546

• Adjustment
  1. Removed the highly collinear variables (Risk?)
  2. Remeasure the collinear variables (e.g., EFA, CFA, IRT)

Problems of using the first adjustment

1. Removing useful information
2. Which variable to remove

2 Heteroscedasticity

2.1 What’s Heteroscedasticity

• Homo/heteroscedasticity ≈ Homo/heterogeneity (of variance)
  • var(ε_i|X) = σ²
  • Violation: \(var(\epsilon_i|X) = \sigma_i^2, \sigma_i^2 \neq \sigma_j^2, \forall i, j.\)

The word “heteroscedasticity” comes from the Greek, and quite literally means data with a different (hetero) dispersion (skedasis).

Regression: Heteroscedasticity

ANOVA: Heterogeneity (of variance)

2.2 Consequence

\[\begin{align} var(\hat\beta|X) =& \frac{\sum(X_i - \bar X)^2}{[\sum(X_i - \bar X)^2]^2}\sigma_i^2.\\ H_0: var(\hat\beta|X) =& \sum(X_i - \bar X)^2;\\ \Leftrightarrow\frac{\sum(X_i - \bar X)^2}{[\sum(X_i - \bar X)^2]^2}\sigma_i^2 =& \frac{\sigma^2}{\sum(X_i - \bar X)^2},\\ \sum(X_i - \bar X)^2\frac{\sigma_i^2}{\sigma^2} =& \sum(X_i - \bar X)^2 \end{align}\]

• Consequence
  • Unbiased, expected consistent, but inefficient
  • Only when \(\frac{\sigma_i^2}{\sigma^2} = 1\), the estimate is efficient.
• Diagnosis: Breusck-Pagan, Goldfeld-Quandt, Park, White, etc.

2.3 Diagnosis: Breusck-Pagan(-Godfreg) test

1. Assuming \(\boldsymbol{\sigma_i^2} = \boldsymbol{X\delta} + \boldsymbol{v_i}\)
2. Regress \(\hat \epsilon_i^2\) on \(X_1, X_2, \cdots, X_{k - 1}\),
   • \(H_0: \sigma_1 = \sigma_2 =\cdots=\sigma_{k - 1} = 0\)
   • \(H_1\): At least one of the σ≠0
   • Statistics \(F_{k - 1, n - k}\)

2.4 Diagnosis: Goldfeld-Quandt test

1. Divide the domain of X to three parts, \(\frac{n -c}{2}, c, \frac{n -c}{2}\).¹
2. Regress Y on X as usual for each part.
3. Calculate \(\sigma_L^2\) from part 1 and \(\sigma_H^2 = \frac{\sum^n_{i = \frac{n + c}{2} + 1}}{\frac{n - c}{2} - k}\hat \epsilon_i^2\) from part 3.
4. Then test \(H_0: \sigma_L^2 = \sigma_H^2\)
   • Statistics: \(\frac{\sigma_L^2}{\sigma_H^2}\sim F_{\frac{n - c - 2k}{2}, \frac{n - c - 2k}{2}}\)

1. Theoretically, the smaller the c is, the more power the test has. However, too small c may also cause the difference between low and high groups difficult to be identified.

2.5 Diagnosis: Park test

Ocular-inspection test: Use the scalar points of \(\hat \epsilon_i^2\) against \(X_i\)¹

Regress \(ln(\hat \epsilon_i^2)\) on some \(X_{ki}\): \[ln(\hat \epsilon_i) = \hat\delta_0 + \hat\delta_1ln(X_{ki}) + \hat\gamma_i.\]

Do the t-test of coefficient in \(ln(X_i)\): \(H_0: \hat\delta_1 = 0.\)

Ocular: 望远镜

1. Or \(\hat \epsilon_i^2\) against \(Y_i\)

2.6 Diagnosis: White Test

1. Regress Y on X to get \(\hat \epsilon_i^2\)
2. Regress \(\hat \epsilon_i^2\) on \(\sum X_i + \sum X_i^2 + \sum X_iX_j\) or \(\hat Y_i\).
3. Test this model against a model with no variables (F test).

2.7 When to use which

• Choose based on data conditions
  • E.g., BP test examines X_i affects the squared residuals, GQ tests differences between subsamples ⇒ Not yield the same result

2.8 Statistics of X1 and X2 from different tests

    Goldfeld-Quandt test

data:  m_fit
GQ = 9.0616, df1 = 4997, df2 = 4997, p-value < 2.2e-16
alternative hypothesis: variance increases from segment 1 to 2

    Goldfeld-Quandt test

data:  m_fit
GQ = 0.94963, df1 = 4997, df2 = 4997, p-value = 0.9661
alternative hypothesis: variance increases from segment 1 to 2

    studentized Breusch-Pagan test

data:  m_fit
BP = 1599.8, df = 1, p-value < 2.2e-16

    studentized Breusch-Pagan test

data:  m_fit
BP = 1.9111, df = 1, p-value = 0.1668

• If you have a specific hypothesis about the variance changes (e.g., increasing with the level of an explanatory variable), the Goldfeld-Quandt or Park test might be appropriate.
• If you’re unsure of the form of heteroscedasticity or have multiple Xs, the Breusch-Pagan or White test would be more suitable.

2.9 Fixed Effect

Least Square with Dummy Variables (LSDV)

• Type I: \(Y_{it} = \beta_1X_{it} + \alpha_i + \epsilon_{it},\) in which α is unit-specific mean differences (unit fixed effect).
• Type II: \(Y_it = \sum^T_{t = 1}\delta_tD_{t_i} + \beta_1X_{it} + \epsilon_{it},\) in which \(\delta_tD_{t_i}\) is the fixed effect for time.
• Issues for using LSDV
  1. Adding a lot of variables (risk of using out of the d.f.);
  2. Additional high multicollinearity;
  3. Losing time invariant variables;
  4. Inefficient estimates of FE on binary and dependent variables.

For the last two concerns, considering using duration models.

2.10 Alternatives corrections

1. Demeaning
2. Fixed effect vector decomposition (FEVD)

Cons:

1. Can’t offer sufficient information (heterogeneous to what extent?).
2. Difficult to calculate substantive effects, e.g., first differences.

⇓

Multilevel modeling (MLM)

2.11 MLM

• Target: Heteroscedasticity can be theoretically meaningful, esp., when researchers focus on how the variance changes.
• Also called as hierarchical linear models (HLM), mixed models, random coefficient models, and variance component models

e.g., Assuming var(ε_i) = σ²e^γ₁X_2i.

If γ₁ = 0, then the model is homoscedastistic; otherwise, one has to use MLE to test H₀: γ₁ = 0.

Modeling variance:

Random Intercept (Two-Level)

\[\begin{align} Y_{ij} = \beta_{0j}& + \beta_{1j}X_{ij} + \epsilon_{ij}, \epsilon_{ij}\sim N(0, \sigma^2)\\ \beta_{0j}& = \gamma_{00} + \gamma_{01}Z_j + u_{0j}, u_{0j}\sim N(0, \tau^2) \end{align}\]

Z is the group indicator.

Intraclass correlation: \(\rho = \frac{\tau^2}{\sigma^2 + \tau^2}.\)

Represent the similarity between individual and group level. 0 means no group effect.

2.12 Random-Slop MLM

\[\begin{align} Y_{ij} = &\beta_{0j} + \beta_{1j}X_{ij} + \epsilon_{ij}\\ &\beta_{0j} = \gamma_{00} + \gamma_{01}Z_j + u_{0j}\\ &\beta_{1j} = \gamma_{10} + u_{1j}.\\ \text{Assume}\left( \begin{array}{c} u_{0j}\\ u_{1j}\end{array}\right)&\sim BVN\left[\left(\begin{array}{c} 0\\0\end{array}\right), \left(\begin{array}{cc} \tau_0^2 & \tau_0\tau_1\\ \tau_0\tau_1 & \tau_1^2 \end{array}\right)\right] \end{align}\]

BVN: Bivariate Normal

The two variables in a bivariate normal are both are normally distributed, and they have a normal distribution when both are added together. Visually, the bivariate normal distribution is a three-dimensional bell curve.

1. Identity: τ₀² = τ₁² = 1;
2. Exchangeable: τ₀²= τ₁²;
3. Independent: τ₀²τ² = 0;
4. Unstructured: no specific relationships is assumed for the τs.

Covariate Matrix

2.13 Alternative Methods

• White covariance matrix
• Weighted least square
• Feasible Generalized Linear Squares (FGLS)
• “Sandwich” Matrix
• Robust Standard Error

2.14 Take-home point

Appendix

2.15 White covariance matrix

(a.k.a., Heteroscedasitic-consistent covariance matrix).

\[E(\epsilon_i|X) = E\{\sqrt{[\epsilon_i - E(\epsilon_i)]^2}|X\}\]

Since E(ε_i|X) = 0 by assumption, we can estimate \(\epsilon_i\) with \(\hat \epsilon_i^2\) and estimate \(var(\hat\beta_1|X)\) with \(\frac{\sum(X_i - \bar X)^2}{[\sum(X_i - \bar X)^2]^2}\hat \epsilon_i^2.\)

The estimates turn out to be biased, but converge asymptotically in n to the true distribution. (Means?)

Means it’s fine when N is large enough.

2.16 Weighted least square

Reduce the substantive effect of the change in X by \(\frac{1}{\sigma_i}\), in order to squeeze the value towards the middle: If σ_i is known,

\[\begin{align} \frac{Y_i}{\sigma_i} =& \frac{\beta_0}{\sigma_i} + \frac{\beta_1}{\sigma_i}X_i + \frac{\epsilon_i}{\sigma_i};\\ \text{Then, } var(\frac{\epsilon_i}{\sigma_i}) =& \frac{1}{\sigma_i^2}var(\epsilon_i) = \frac{\sigma_i^2}{\sigma_i^2} = 1;\\ \Rightarrow Y_i^* =& \beta_0X^*_{0i} + \beta_1X^*_{1i} + \epsilon_i^*. \end{align}\]

The last equation is homoscedastistic. However, in most cases, we don’t know σ_i ⇒ var(ε_i)∼ X_1i, i.e., var(ε_i) = σ_i² = σ²X_i = h_iσ².

2.17 Weighted least square (continued)

Given the goal var(ε^*_i) = var(ε^*_j), ∀ i,j.

\[\begin{align} \epsilon^*_i =& \frac{\epsilon_i}{\sqrt{h_i}},\\ var(\frac{\epsilon_i}{\sqrt{h_i}}) =& \frac{var(\epsilon_i)}{h_i} = \frac{h_i\sigma^2}{h_i} = \sigma^2,\\ \Rightarrow Y_i^* =& \frac{Y_i}{\sqrt{h_i}}; X_{0i}^* = \frac{1}{\sqrt{h_i}}; X_{1i}^* = \frac{X_{1i}}{\sqrt{h_i}}, \text{ assuming } X_i\in R^+. \end{align}\]

Tip

In practice, there are different ways to estimate the weight, one need to carefully choose the proper one.

2.18 Feasible Generalized Linear Squares (FGLS)

\[\begin{align} \boldsymbol{Y} =& \boldsymbol{X\beta} + \boldsymbol{\epsilon}; \\ var(\boldsymbol{\epsilon}) =& \Omega_{n\times n} = \left( \begin{array}{cccc} \sigma_1^2 & 0 & \cdots & 0\\ 0 & \sigma_2^2 & \cdots & 0\\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & \sigma_n^2\\ \end{array}\right),\\ \text{Then, } \boldsymbol{\hat\beta_{GLS}} =& (\boldsymbol{X'\Omega X})^{-1}(\boldsymbol{X'\Omega Y}).\\ \text{Let }\boldsymbol{H}: \boldsymbol{\Omega} = \boldsymbol{HH^{-1}},\text{then, } \boldsymbol{H^{-1}Y} =& \boldsymbol{H^{-1}X\beta} + \boldsymbol{H^{-1}\epsilon}, \boldsymbol{H^{-1}\epsilon} = \boldsymbol{H^{-1}(H^{-1})'};\\ var(\boldsymbol{\epsilon})=& (\boldsymbol{HH'})^{-1}\boldsymbol{\Omega} = \boldsymbol{\Omega}^{-1}\boldsymbol{\Omega} = \boldsymbol{I},\\ var(\boldsymbol{\hat\beta_{GLS}}) =& (\boldsymbol{X'X})^{-1}(\boldsymbol{X'\Omega X})(\boldsymbol{X'X})^{-1}. \end{align}\]

• If there’s no heteroscedasticity, \(\boldsymbol{\Omega} = \sigma\boldsymbol{I}.\)
• If there is heteroscedasticity and \(\boldsymbol{\Omega}\) is known, then WLS (a special type of GLS).
• If \(\boldsymbol{\Omega}\) is unknown, run \(\boldsymbol{Y} = \boldsymbol{X\beta} + \boldsymbol{\epsilon}\) to estimate \(\boldsymbol{\Omega}\) with \(\boldsymbol{\hat\Omega} = \boldsymbol{\hat\epsilon\hat\epsilon'}\)

Warning

NB: This method does not get SE, and also biased for small N.

2.19 “Sandwich” Matrix

FGLS is a type of “sandwich” estimator.

In a more general view, let \(\boldsymbol{Q} = \boldsymbol{X'X}\) and

\[\begin{align} \boldsymbol{Q} = \left( \begin{array}{cccc} \sigma_1^2 & 0 & \cdots & 0\\ 0 & \sigma_2^2 & \cdots & 0\\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & \sigma_n^2\\ \end{array}\right) \end{align}\]

Then for regular OLS, \(var(\beta) = \sigma^2(\boldsymbol{X'X})^{-1} = \sigma^2\boldsymbol{Q}^{-1}\).

But when heteroscedasticity occurs, \(var(\boldsymbol{\beta}|X)\neq \sigma^2\boldsymbol{Q}^{-1}\).

Instead, let \(\boldsymbol{G} = \boldsymbol{X'GX}\), then \(var(\boldsymbol{\beta}|X) = \boldsymbol{Q^{-1}GQ}^{-1}.\)

2.20 Heteroscedasticity in TSCS Data

• TSCS: Time-series cross-sectional data.
  • Sometimes, they are called “pooled data”, “panel data”(WRONG!)
• Reason of heterosceasticity
  • Spatical heterogeneity
  • Temporal heterogeneity
• Deal with heteroscedasticity in TSCS data:
  1. Robust standard error
  2. Fixed effect
  3. Multilevel modeling

2.21 Robust Standard Error

• Sandwich SE
  • In the FGLS version, the “meat” can be identified only when T > N, and before this, the autocorrelation (\(cov(\epsilon_i, u_j|X_i, X_j) = 0, \forall i, j\)) has to be eliminated.
• Panel-corrected SE
  • An alternative version of the sandwich SE. Ω is clustered by time periods: \(\hat\Omega_{ij} = \frac{\sum^T_{t = 1}\epsilon_{it}\epsilon_{jt}}{T}\)
• Cluster SE
  • Adjust SE to account for correlations within clusters.
• Within-between model
  1. Run separate models in each group
  2. Run an aggregate regression

Let \(\bar Y_i = \frac{\sum^n_{i=1}Y_{it}}{n_i}, \bar X_i = \frac{\sum^n_{i=1}X_{it}}{n_i}, \bar \epsilon_i = \frac{\sum^n_{i=1}\epsilon_{it}}{n_i}\). Then,

\[\begin{align} \bar{Y_i} =& \beta_0 + \beta_1\bar{X_i} + \epsilon_i + \epsilon_i\\ Y_{it} - \bar{Y_i} =& (\beta_0 - \beta_0) + \beta_1(X_{it} - \bar{X_i}) + (\epsilon_{it} - \epsilon_i)\\ \hat{Y_i} =& \beta_1\hat{X_{it}} + \hat{\epsilon_{it}}, \text{(a.k.a., the within model)}\\ Y_{it} =& \beta_0 + \beta_1(X_{it} - \bar{X_i}) + \beta_2\bar{X_i} + \epsilon_{it}, \text{(a.k.a., the between model)} \end{align}\]

Stretch

https://www.bilibili.com/video/BV1E54y1r76o/?buvid=Y045340D2FD9FA3A46419EEFE4578279ECBD&from_spmid=main.space-search.0.0&is_story_h5=false&mid=7RnjBONLRMus4FQZFBWD2g%3D%3D&p=2&plat_id=114&share_from=ugc&share_medium=iphone&share_plat=ios&share_session_id=79FE08C9-F2B1-4C18-A16C-B31179817B42&share_source=WEIXIN&share_tag=s_i&timestamp=1728732278&unique_k=skSWw7j&up_id=390316092&vd_source=f38aeefd0d38cecba9017eeee43e71c8

Meditation

松茸的世界：5分钟正念冥想-自信之心

https://www.xiaoyuzhoufm.com/episode/65b752ed0bef6c207457f51b?s=eyJ1IjogIjYwNDA1OGJiZTBmNWU3MjNiYmY3MjZiZSJ9