The R Project SVN R-packages

library(nlme)

## PR#17712: nlme() with a fixed variance function
set.seed(42)
d <- data.frame(
    obs = rnorm(137), # prime, so recycling would fail
    group = factor(rep_len(c("A", "B"), 137))
)
fit_wt <- function (weights = NULL)
    nlme(obs ~ b, fixed = b ~ 1, random = b ~ 1 | group, start = c(b = 0),
         data = d, weights = weights,
         control = nlmeControl(returnObject = TRUE)) # possibly small PNLS steps

m_unweighted <- fit_wt(NULL)
stopifnot(all.equal(fixef(m_unweighted), c(b = mean(d$obs)), tolerance = 0.05))
## equivalent fit with a constant variance function:
m_varIdent <- fit_wt(varIdent())
## nlme <= 3.1-164 failed with Error in recalc.varFunc(object[[i]], conLin) :
##   dims [product 12] do not match the length of object [13]
stopifnot(all.equal(logLik(m_unweighted), logLik(m_varIdent)))

fit0 <- fit_wt(varExp(form = ~obs))
## equivalent fit with the coefficient **fixed** at its estimate:
vf0 <- fit0$modelStruct$varStruct
fit1 <- fit_wt(varExp(form = ~obs, fixed = unname(coef(vf0))))
## fitting failed in nlme <= 3.1-164 (as above)
stopifnot(all.equal(logLik(fit0), logLik(fit1),
                    check.attributes = FALSE)) # df differs by 1
## equivalent fit using the same variance weights via **varFixed**:
##d$wt <- varWeights(vf0)  # beware of different (internal) ordering
d$wt <- fit0$sigma / attr(fit0$residuals, "std")
fit2 <- fit_wt(varFixed(~1/wt^2))
## fitting failed in nlme <= 3.1-164 (as above)
stopifnot(all.equal(logLik(fit0), logLik(fit2),
                    check.attributes = FALSE)) # df differs by 1