2
$\begingroup$

I am new to both mixed effect and Additive models so I'm sorry if the answer here is trivial.

I have data collected on several metabolic chemicals (M1,M2...), covariates (time,Race,Gender...) and disease state (D,D.binary). I'm trying to generate a GAMM based on variables selected from a GEE variable selection.

Data:

  • 8 cases, 51 matched controls
  • approximately 10 time points from each subject
  • ~630 observations
  • M1,M2...M3 are metabolites many of which are formed from common parts, Metabolite levels are correlated in that they are competing for the same component parts
  • Covariates stratify the subjects into subgroups

Here is my model as it is now:

> b = gamm(D.binary ~ Time  + s(M1)  , 
      random = list(ParticipantID = ~ 1 + Time),  niterPQL=50,
      data = NEC_data, family=binomial(link="logit")) 

 Maximum number of PQL iterations:  50 
 iteration 1
 iteration 2
 ...
 iteration 49
 iteration 50
 Warning message:
 In gammPQL(y ~ X - 1, random = rand, data = strip.offset(mf), family = family,  :
  gamm not converged, try increasing niterPQL

> plot(b$gam,pages=1)

enter image description here

> summary(b$lme) # details of underlying lme fit
Linear mixed-effects model fit by maximum likelihood
 Data: data 
   AIC  BIC logLik
  -160 -124     88

Random effects:
 Formula: ~Xr - 1 | g
 Structure: pdIdnot
          Xr1   Xr2   Xr3   Xr4   Xr5   Xr6   Xr7   Xr8
StdDev: 0.812 0.812 0.812 0.812 0.812 0.812 0.812 0.812

 Formula: ~1 + Time | ParticipantID %in% g
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev  Corr  
(Intercept) 5.68324 (Intr)
Time        0.50739 -0.92 
Residual    0.00691       

Variance function:
 Structure: fixed weights
 Formula: ~invwt 
Fixed effects: list(fixed) 
                   Value Std.Error  DF t-value p-value
X(Intercept)       -2.81     0.729 573   -3.86  0.0001
XTime              -0.15     0.065 573   -2.30  0.0220
Xs(M1)Fx1 -1.60     0.066 573  -24.29  0.0000
 Correlation: 
                   X(Int) XTime  
XTime              -0.920       
Xs(M1)Fx1  0.004  0.000

Standardized Within-Group Residuals:
    Min      Q1     Med      Q3     Max 
-2.3472 -0.0692 -0.0117  0.0305 20.7271 

Number of Observations: 636
Number of Groups: 
                   g ParticipantID %in% g 
                   1                   61 
> summary(b$gam) # gam style summary of fitted model

Family: binomial 
Link function: logit 

Formula:
NEC.binary ~ Time + s(M1)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -2.8135     0.7289   -3.86  0.00013 ***
Time         -0.1495     0.0651   -2.30  0.02188 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
               edf Ref.df     F p-value    
s(M1) 4.1    4.1 14913  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.0872   
  Scale est. = 4.7815e-05  n = 636
> anova(b$gam) 

Family: binomial 
Link function: logit 

Formula:
NEC.binary ~ Time + s(M1)

Parametric Terms:
    df    F p-value
Time  1 5.28   0.022

Approximate significance of smooth terms:
               edf Ref.df     F p-value
s(M1) 4.1    4.1 14913  <2e-16
> gam.check(b$gam)

enter image description here

I suspect I may have messed up something fairly basic since M1 is the most obvious discriminator of the disease state. It is significant (as it should be) but the correlation is very low. Also, obviously, the model didn't converge (even when I increased iterations from 20->50). And finally the check plots look pretty outrageous

Question

Have I made a basic syntax error? Is there some malicious component in my model I'm over looking? Any help would be greatly appreciated.

Further work

I would like to add another metabolite (M2) to the model and 2 covariates (Birthweight and Race). When I add M2 to the model I get an non-convergence error:

> b = gamm(D.binary ~ Time  + s(M1) + s(M2) , 
      random = list(ParticipantID = ~ 1 + Time),  niterPQL=20, correlation = corLin(),
      data = NEC_data, family=binomial(link="logit"))

 Maximum number of PQL iterations:  20 
iteration 1
iteration 2
Error in lme.formula(fixed = fixed, random = random, data = data, correlation = correlation,  : 
  nlminb problem, convergence error code = 1
  message = false convergence (8)

Any advice about moving into the multidimensional space would also be appreciated.

Addition

I also tried this model with the discrete disease classification (control: 0,1 disease: 2,3) and poisson noise.

> b = gamm(NEC ~ DPP  + s(DSLNT_ug.mL)  , 
+      random = list(ParticipantID = ~ 1 + DPP),  niterPQL=20,
+      data = NEC_data, family=poisson) 

 Maximum number of PQL iterations:  20 
iteration 1
iteration 2
...
iteration 19
iteration 20
Error in solve.default(pdMatrix(a, factor = TRUE)) : 
  system is computationally singular: reciprocal condition number = 3.13906e-19
Improve this question
$\endgroup$
10
  • $\begingroup$ It's hard to know if this is what's happening at a distance, but I ran into a very similar problem awhile ago. In particular, suppose you have subjects whose outcomes look like this: 0, 0, 0, 1, 1, 1, ...Then if you fit the individual's intercept and slope, you will get $\hat \beta_0 = -\infty$, and $\hat \beta_1 = \infty$. A mixed effects model may help "pull" these estimates in, but it's not hard to imagine the problems this creates if this appears frequently in your dataset. $\endgroup$
    – Cliff AB
    Commented Jun 4, 2015 at 18:25
  • $\begingroup$ In my case, I switched to a GEE model; on the population level, I did not see infinite slopes, and I was really more interested in population level effects anyways. $\endgroup$
    – Cliff AB
    Commented Jun 4, 2015 at 18:25
  • 1
    $\begingroup$ @CliffAB Our collaborating biostatistician advised against using GEE for a predictive model. The explanation was lost in translation but that is why I'm attempting to get this answer from GAMM. I get a different error when I run it with the discrete disease classification (edited into the question) indicative on an un-invertible design matrix due to too many highly correlated variables (how can that be, I only have one fixed effect?). $\endgroup$
    – user3030872
    Commented Jun 4, 2015 at 19:36
  • $\begingroup$ Ah yes, if you want a predictive model, GEE's are not the way to go. Different study goals. I should have realized that given that you are using GAMM's. But while I don't have an answer for you, it might be worth reviewing your data to see if this trend of potential infinite slopes exists. If so, you might want to think about parameterizing the problem in a different way. $\endgroup$
    – Cliff AB
    Commented Jun 4, 2015 at 19:45
  • $\begingroup$ @CliffAB I don't see the infinite slopes in my summaries here. When you say "reviewing your data to see if this trend of potential infinite slopes exists." Do you mean remake these models with each of my variables and see if the problem continues? Also, when you say "If so, you might want to think about parameterizing the problem in a different way." Can you expand on what parameters I should be playing with? thanks. $\endgroup$
    – user3030872
    Commented Jun 4, 2015 at 23:24

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.