[R-sig-ME] mgcv gamm problem with factor-smooth interaction and AR1 errors

[Phillip Alday]

Without really trying to understand your models and whether my
suggestion does the same thing, I found that this works:

m6 = gamm(x~s(index, by=subject, bs='fs'),
correlation=corAR1(form=~index|subject), data=df, method = 'REML')

I just looked at the differences between lines and saw that "subject"
was specified as a named argument in some lines, but not others.

Phillip

On 04/23/2018 05:46 PM, Mathew Guilfoyle wrote:
> 
> I've posted this question to the main R help list without any joy - I wonder if someone here might be able to help...
> 
> I'm trying to fit a mgcv::gamm model including random smooths and an autocorrelation term but am getting a consistent error.  I get the same issue with my real data and the toy data in the example below.  As far as I can see the model is specified correctly and the mgcv man/help pages imply this sort of model should be possible.
> 
> Any ideas?
> 
> Many thanks
> Mathew
> 
> library('mgcv')
> set.seed(1)
> #generate some data
> df = data.frame(index=rep(1:10,5), x=runif(50,0,1), subject = as.factor(sort(rep(1:5,10))))
> 
> ##models M1-M4 all fit without issue, but M5 fails
> 
> # random intercept
> m1 = gamm(x~s(index), random=list(subject=~1), data=df, method = 'REML')
> 
> #factor interaction, random intercept, AR errors
> m2 = gamm(x~s(index, by=subject), random=list(subject=~1), correlation=corAR1(form=~index|subject), data=df, method = 'REML')
> 
> #factor interaction, random intercept and slope, AR errors
> m3 = gamm(x~s(index, by=subject), random=list(subject=~index), correlation=corAR1(form=~index|subject), data=df, method = 'REML')
> 
> # the 'fs' smooth on its own works ok
> m4 = gamm(x~s(index, subject, bs='fs'), data=df, method = 'REML')
> 
> #combination of 'fs' smooth and AR errors generates the error: "Error in matrix(0, size.cg[i], size.cg[i]) : object 'size.cg' not found"
> m5 = gamm(x~s(index, subject, bs='fs'), correlation=corAR1(form=~index|subject), data=df, method = 'REML')
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>