variance function on a mean polynomial (nlme R package)

Question

I want to fit a variance function that involves polynomial terms of the mean

I know that I can fit the variance function

using varExp() (without arguments). Reading the documentation of varExp() it defaults to ~ fitted(.), however using varExp(~ fitted(.)) gives an error, therefore varExp(~ fitted(.)^2) is not working.

Here is a code example in which I am working:

library(pacman)
p_load(tidyverse)
p_load(MASS)
p_load(nlme)
set.seed(123)
n<-1000

z0<-data.frame(days=runif(n,-3,3))

z0_mat <- model.matrix(~ days, z0)
betas<-c(10,.85)
a0<-1.5
b0<-8
b1<- -0.2
x_target<-60
alpha<-(log(b0)-log(a0))/x_target
z0$weight <- z0_mat%*%betas+rnorm(n,0,a0*exp(alpha*z0_mat[,2]+b1*z0_mat[,2]^2))
z0 %>% ggplot(aes(x=days,y=weight))+geom_point()

m01<-gls(weight~days,weights=varComb(varExp(form=~days),varExp(form=~days^2)),data=z0)
intervals(m01,which = "var-cov")
#> Approximate 95% confidence intervals
#> 
#>  Variance function:
#>                lower        est.       upper
#> A.expon -0.007686726  0.01778161  0.04324994
#> B.expon -0.212892581 -0.19686001 -0.18082744
#> 
#>  Residual standard error:
#>    lower     est.    upper 
#> 1.394629 1.487835 1.587270

^{Created on 2022-09-13 by the reprex package (v2.0.1)}

in m01 I implemented a variance function on the covariate days,

but I would want to implement it on the mean. In this example this approach would be enough, however with more covariates in (the mean structure of) the model it would be better to use the mean rather than covariates.

Answer 1

This runs without throwing an error

v1  <- varComb(
    varExp(form=~fitted(.)),
    varExp(form=~I(fitted(.)^2)))
m01 <- gls(weight~days, weights= v1, data = z0)

• It is always confusing that form is not the first argument to the var* functions, thus varExp(fitted(.)) doesn't do what you think it does ...
• I thought you would need the I(...) in the quadratic term, but apparently you don't ...
• However, the reported residual std error is vanishingly small — I'm not sure what's up with that.