glsD {Design} | R Documentation

## Fit Linear Model Using Generalized Least Squares

### Description

This function fits a linear model using generalized least squares. The errors are allowed to be correlated and/or have unequal variances. `glsD` is a slightly enhanced version of the Pinheiro and Bates `glsD` function in the `nlme` package to make it easy to use with the Design library 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).

### Usage

    glsD(model, data, correlation, weights, subset, method, na.action,
        control, verbose, B=0, dupCluster=FALSE, pr=FALSE,
        opmeth=c('optimize','optim'))

    ## S3 method for class 'glsD':
    print(x, digits=4, ...)

### Arguments

`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.fail`) causes `gls` to print an error message and terminate if there are any incomplete observations.
`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 `glsD` 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
`opmeth` | specifies whether the `optimize` or the `optim` function is to be used for optimization
`x` | the result of `glsD`
`digits` | number of significant digits to print
`...` | ignored

### Value

an object of classes `glsD`, `Design`, and `gls` representing the linear model fit. Generic functions such as `print`, `plot`, 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. `glsD` returns the following components not returned by `gls`: `Design`, `assign`, `formula`, `opmeth` (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).

### Author(s)

Jose Pinheiro jcp@research.bell-labs.com, Douglas Bates bates@stat.wisc.edu, Frank Harrell f.harrell@vanderbilt.edu, Patrick Aboyoun aboyoun@insightful.com

### References

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

### See Also

`gls` `glsControl`, `glsObject`, `varFunc`, `corClasses`, `varClasses`

### Examples

    ## 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')
    # library(nlme)  # S-Plus: library(nlme3) or later
    f <- glsD(y ~ pol(tim,2) + trt, correlation=corCAR1(form= ~tim | id),
              data=d, B=B)
    f
    f$var      # bootstrap variances
    f$varBeta  # original variances
    summary(f)
    anova(f)
    plot(f, tim=NA, trt=NA)
    # v <- Variogram(f, form=~tim|id, data=d)
    ## End(Not run)

[Package Contents]