拟合GLM(inverse.gaussian)到AR(1)模拟数据的自回归估计问题求助

Hey there, let's break down this tricky inverse Gaussian GLM + AR(1) problem you're tackling—dealing with skewed groups and autoregressive effects can get messy, so let's start with the key spots to investigate:

1. Align Your Theoretical Skew-Variance Relationship with the Inverse Gaussian Distribution #

First, let's ground ourselves in the inverse Gaussian (IG) distribution's inherent properties:

• Skewness: (3\sqrt{\mu/\lambda}) (where (\mu) = mean, (\lambda) = precision parameter)
• Variance: (\mu^3/\lambda)

Your theory states higher skew should correspond to smaller variance, but the IG's variance scales with (\mu^3/\lambda). For this relationship to hold, you need to adjust (\mu) and (\lambda) in tandem: if you increase skew (by raising the ratio (\mu/\lambda)), you have to decrease (\mu) enough to offset the (\mu^3) term in the variance formula. For example:

• Group 1 (low skew): (\mu=3), (\lambda=9) → skew = (3\sqrt{3/9}=1.73), variance = (27/9=3)
• Group 4 (high skew): (\mu=1), (\lambda=1) → skew = (3\sqrt{1/1}=3), variance = (1/1=1)

If your simulated data doesn't follow this parameterization, your model results will contradict your theory right out the gate. Double-check how you're generating each group's data—use tools like moments::skewness() and var() in R to verify sample stats match your theoretical targets.

2. Don't Ignore AR(1) Autocorrelation in Standard GLMs #

Standard GLMs assume independent observations, but your AR(1) data has built-in autocorrelation between consecutive measurements. Ignoring this will:

• Bias standard errors (leading to incorrect p-values)
• Make parameter estimates inefficient
• Possibly skew your autoregressive effect estimates entirely

Instead, use methods designed for correlated data:

• Generalized Estimating Equations (GEE): Specify an AR(1) working correlation structure. Great for population-averaged effects and robust to minor correlation misspecification. In R, use geepack::geeglm().
• Generalized Linear Mixed Models (GLMM): Add a structured covariance matrix (like AR(1)) to account for within-subject autocorrelation. Packages like glmmTMB or nlme support this for IG models.
• Bayesian Models: Use brms or rstanarm to explicitly model the AR(1) structure alongside the IG likelihood—this gives you flexible uncertainty estimates and handles complex relationships easily.

3. Choose the Right Link Function for Autoregressive Effect Interpretation #

The IG GLM's link function dictates how you interpret the autoregressive (X_{t-1} → X_t) effect:

• Inverse link (default in R's glm()): (1/\mu_t = \beta_0 + \beta_1 X_{t-1}) → A unit increase in X_{t-1} decreases the mean of X_t by (\beta_1) times the inverse of the current mean.
• Log link: (\log(\mu_t) = \beta_0 + \beta_1 X_{t-1}) → A unit increase in X_{t-1} multiplies the mean of X_t by (e^{\beta_1}) (multiplicative effect).
• Identity link: (\mu_t = \beta_0 + \beta_1 X_{t-1}) → Linear additive effect, but risky since IG means are positive and the link doesn't enforce this.

Pick a link that matches your theoretical expectations of how X's past values influence its current values. Also, include group-specific intercepts (and maybe slopes) to separate skew-related mean/variance differences from the core autoregressive effect.

4. Diagnose Model Fit for Hidden Issues #

Once you've adjusted your model to handle autocorrelation, run diagnostics to check for problems:

• Plot residual autocorrelation (ACF plots) to confirm the AR(1) structure is being captured—if residuals still show strong autocorrelation, your correlation model is misspecified.
• Check for convergence warnings (common with GLMMs/Bayesian models). If the model won't converge, try starting values closer to plausible estimates or simplify the model temporarily.
• Compare parameter estimates across groups: Do high-skew groups show the smaller autoregressive effects you expect? If not, revisit your simulation setup or theoretical assumptions about how skew and shared variance interact.

If you can share specific code snippets, model output, or the exact issues you're hitting (e.g., convergence failures, unexpected coefficient signs), we can dive even deeper!

内容的提问来源于stack exchange，提问作者Jeroen Mulder