[R-sig-ME] lme() - Error in recalc.corAR1

[Santosh Srinivas]

Thank you very much, Dr. Bolker for such a detailed and quick response. Very much appreciate the troubleshooting steps

I can confirm that in our earlier iteration when x and y were grand-mean entered (just as you have suggested below), the  lme model converged without any issues.

However, we tried to fit the same model with this data where x and y are now centered within user_id. This change was done following the recommendations in the below paper:
Enders, Craig K., and Davood Tofighi. 2007. �Centering Predictor Variables in Cross-Sectional Multilevel Models: A New Look at an Old Issue.� Psychological Methods 12(2): 121�38. Doi:10.1037/1082-989X.12.2.121. https://osf.io/mt5z3/download

Not sure why a lme model with centering at grand mean (CGM) works, but the one with centering within cluster (CWC) fails.

We are currently trying glmmTMB. It has been running for last 2-3 hours with CWC data. Hopefully the model converges.

Many thanks, again!

Regards
Santosh
________________________________
From: R-sig-mixed-models <r-sig-mixed-models-bounces_at_r-project.org_santoshbsrinivas using duck.com> on behalf of Ben Bolker <bbolker_at_gmail.com_santoshbsrinivas using duck.com>
Sent: Thursday, January 2, 2025 5:23 PM
To: r-sig-mixed-models using r-project.org <r-sig-mixed-models_at_r-project.org_santoshbsrinivas using duck.com>
Subject: Re: [R-sig-ME] lme() - Error in recalc.corAR1

DuckDuckGo was unable to verify sender identity

   Incomplete solution: scaling x and y seems to help (this is a
surprisingly common solution for numerical problems, which lme is
particularly prone to for large  data sets.

   However, I didn't try to fit the full model (500K observations); I
used a random subsample of 10,000 points. You should try gradually
increasing the size of the subsample to see how the problem scales
(i.e., will you have to wait 3 minutes or 3 days for the full answer?)
and to see if you run into any other problems

library(dplyr) ## data manipulation -- could also be done in base R
library(nlme)
library(glmmTMB)
library(ggplot2); theme_set(theme_bw())

set.seed(101)

dd <- read.csv("r-sig-mixed-02jan.csv")

## add per-user time step indicator, and dummy group (for glmmTMB);
##  scale x and y
dd2 <- (dd
     |> group_by(user_id)
     |> mutate(step = factor(seq(n())))
     |> ungroup()
     |> mutate(xs = drop(scale(x)), ys = drop(scale(y)),
               dummy = factor(1))
)

## subsample
dds <- dd2[sample(nrow(dd2), size = 1e4),]


ggplot(dd, aes(x, y, group = user_id)) + geom_line(alpha = 0.2)


system.time(m <- lme(
     ys ~ xs,
     random = ~xs|user_id,
     correlation = corAR1(form = ~1|user_id),
     weights = varExp(form = ~xs),
     method = 'REML',
     na.action = na.omit,
     data = dds,
     control = lmeControl(
         maxIter = 1000, msMaxIter = 1000,
         msVerbose = TRUE,
         returnObject = TRUE  ## return result even w/ optimizer warning
      )
     )
     )
## 3 seconds

The "singular convergence" warning from nlminb is inscrutable but
possibly ignorable.
https://stackoverflow.com/questions/79110546/glmmtmb-convergence-messages


Looking at the output, it seems that you will probably have to simplify
your model in any case -- there is a correlation of 0.999 between the
random intercepts and slopes. The estimated AR1 coefficient is also very
small


## Try full data set. If the problem scales *linearly* (optimistic!)
## then we should expect this to take ~ 150 seconds ...

system.time(m2 <- update(m, data = dd2))


   glmmTMB works too, and gets similar answers, but is much slower in
this case (6 minutes vs 3 seconds); turning on some parallelization
*might* speed it up.

system.time(m3 <- glmmTMB(ys ~ xs + (xs|user_id) + ar1(0 + step |user_id),
         dispformula = ~0 + xs,  ## don't want a nugget in corAR1
         REML = TRUE,
         data = dds,
         verbose = TRUE)
         )

fixef(m)
fixef(m3)$cond

VarCorr(m)
VarCorr(m3)


On 2025-01-02 9:08 a.m., Santosh Srinivas wrote:
> Dear List Recipients, Hope this email finds you well.
>
> I am writing to request your help the below issue:
> Error in recalc.corAR1(object[[i]], conLin): NA/NaN/Inf in foreign function call (arg 1)
>
> This occurs when I run the following regression model:
>
> m =
>    lme(
>      y ~ x,
>      random = ~x|user_id,
>      correlation = corAR1(form = ~1|user_id),
>      weights = varExp(form = ~x),
>      method = 'REML',
>      na.action = na.omit,
>      data = df)
>
> I have uploaded the data here for reproducibility:
> https://drive.google.com/file/d/1kA5-T4zZ-f33gJgAsJ9AWMlK-YXAm65m/view?usp=sharing
>
> Request your help with how to understand and fix this ussue.
>
> Thanks & regards
> Santosh
>
>
>
>
>
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

--
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Director, School of Computational Science and Engineering
 > E-mail is sent at my convenience; I don't expect replies outside of
working hours.

_______________________________________________
R-sig-mixed-models using r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

	[[alternative HTML version deleted]]