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

Fri Jan 3 03:49:01 CET 2025

   OK. (It would have been helpful to have that context ...)

   A couple of quick thoughts --

* I would definitely suggest some benchmarking on smaller subsets, if 
you have the time.

* I can see why you don't have flexibility in centering, but *scaling* 
(which only changes the scale of the parameters, not anything about the 
conclusions) might also be numerically helpful. (It's unfortunate that 
R's scale() function conflates scaling and centering in its default 
options ...)

* The model fit I got on the smaller data set definitely suggested some 
model simplification would be reasonable ...

   cheers
    Ben Bolker

On 2025-01-02 11:41 a.m., Santosh Srinivas wrote:
> 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]]
> 
> 
> _______________________________________________
> 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.