This function fits a linear model using generalized least
squares. The errors are allowed to be correlated and/or have unequal
variances. `Gls`

is a slightly enhanced version of the
Pinheiro and Bates `gls`

function in the `nlme`

package to
make it easy to use with the rms package and to implement cluster
bootstrapping (primarily for nonparametric estimates of the
variance-covariance matrix of the parameter estimates and for
nonparametric confidence limits of correlation parameters).

For the `print`

method, format of output is controlled by the
user previously running `options(prType="lang")`

where
`lang`

is `"plain"`

(the default), `"latex"`

, or
`"html"`

.

```
Gls(model, data, correlation, weights, subset, method, na.action=na.omit,
control, verbose, B=0, dupCluster=FALSE, pr=FALSE, x=FALSE)
```# S3 method for Gls
print(x, digits=4, coefs=TRUE, title, …)

model

a two-sided linear formula object describing the
model, with the response on the left of a `~`

operator and the
terms, separated by `+`

operators, on the right.

data

an optional data frame containing the variables named in
`model`

, `correlation`

, `weights`

, and
`subset`

. By default the variables are taken from the
environment from which `gls`

is called.

correlation

an optional `corStruct`

object describing the
within-group correlation structure. See the documentation of
`corClasses`

for a description of the available `corStruct`

classes. If a grouping variable is to be used, it must be specified in
the `form`

argument to the `corStruct`

constructor. Defaults to `NULL`

, corresponding to uncorrelated
errors.

weights

an optional `varFunc`

object or one-sided formula
describing the within-group heteroscedasticity structure. If given as
a formula, it is used as the argument to `varFixed`

,
corresponding to fixed variance weights. See the documentation on
`varClasses`

for a description of the available `varFunc`

classes. Defaults to `NULL`

, corresponding to homoscesdatic
errors.

subset

an optional expression indicating which subset of the rows of
`data`

should be used in the fit. This can be a logical
vector, or a numeric vector indicating which observation numbers are
to be included, or a character vector of the row names to be
included. All observations are included by default.

method

a character string. If `"REML"`

the model is fit by
maximizing the restricted log-likelihood. If `"ML"`

the
log-likelihood is maximized. Defaults to `"REML"`

.

na.action

a function that indicates what should happen when the
data contain `NA`

s. The default action (`na.omit`

) results
in deletion of observations having any of the variables of interest missing.

control

a list of control values for the estimation algorithm to
replace the default values returned by the function `glsControl`

.
Defaults to an empty list.

verbose

an optional logical value. If `TRUE`

information on
the evolution of the iterative algorithm is printed. Default is
`FALSE`

.

B

number of bootstrap resamples to fit and store, default is none

dupCluster

set to `TRUE`

to have `Gls`

when
bootstrapping to consider multiply-sampled clusters as if they were
one large cluster when fitting using the `gls`

algorithm

pr

set to `TRUE`

to show progress of bootstrap resampling

x

for `Gls`

set to `TRUE`

to store the design matrix
in the fit object; otherwise the result of `Gls`

digits

number of significant digits to print

coefs

specify `coefs=FALSE`

to suppress printing the table
of model coefficients, standard errors, etc. Specify `coefs=n`

to print only the first `n`

regression coefficients in the
model.

title

a character string title to be passed to `prModFit`

…

ignored

an object of classes `Gls`

, `rms`

, and `gls`

representing the linear model
fit. Generic functions such as `print`

, `plot`

,
`ggplot`

, and `summary`

have methods to show the results of
the fit. See
`glsObject`

for the components of the fit. The functions
`resid`

, `coef`

, and `fitted`

can be used to extract
some of its components. `Gls`

returns the following components
not returned by `gls`

: `Design`

, `assign`

,
`formula`

(see arguments), `B`

(see
arguments), `bootCoef`

(matrix of `B`

bootstrapped
coefficients), `boot.Corr`

(vector of bootstrapped correlation
parameters), `Nboot`

(vector of total sample size used in each
bootstrap (may vary if have unbalanced clusters), and `var`

(sample variance-covariance matrix of bootstrapped coefficients). The
\(g\)-index is also stored in the returned object under the name
`"g"`

.

The `na.delete`

function will not work with
`Gls`

due to some nuance in the `model.frame.default`

function. This probably relates to `na.delete`

storing extra
information in the `"na.action"`

attribute of the returned data
frame.

Pinheiro J, Bates D (2000): Mixed effects models in S and S-Plus. New York: Springer-Verlag.

`gls`

`glsControl`

, `glsObject`

,
`varFunc`

, `corClasses`

,
`varClasses`

, `GiniMd`

,
`prModFit`

, `logLik.Gls`

# NOT RUN { ns <- 20 # no. subjects nt <- 10 # no. time points/subject B <- 10 # no. bootstrap resamples # usually do 100 for variances, 1000 for nonparametric CLs rho <- .5 # AR(1) correlation parameter V <- matrix(0, nrow=nt, ncol=nt) V <- rho^abs(row(V)-col(V)) # per-subject correlation/covariance matrix d <- expand.grid(tim=1:nt, id=1:ns) d$trt <- factor(ifelse(d$id <= ns/2, 'a', 'b')) true.beta <- c(Intercept=0,tim=.1,'tim^2'=0,'trt=b'=1) d$ey <- true.beta['Intercept'] + true.beta['tim']*d$tim + true.beta['tim^2']*(d$tim^2) + true.beta['trt=b']*(d$trt=='b') set.seed(13) library(MASS) # needed for mvrnorm d$y <- d$ey + as.vector(t(mvrnorm(n=ns, mu=rep(0,nt), Sigma=V))) dd <- datadist(d); options(datadist='dd') f <- Gls(y ~ pol(tim,2) + trt, correlation=corCAR1(form= ~tim | id), data=d, B=B) f AIC(f) f$var # bootstrap variances f$varBeta # original variances summary(f) anova(f) ggplot(Predict(f, tim, trt)) # v <- Variogram(f, form=~tim|id, data=d) nlme:::summary.gls(f)$tTable # print matrix of estimates etc. options(datadist=NULL) # }