getG: Extract R, G or V matrix in a mixed or GLS model

In gmonette/gnew: New functions for spida2

Extract R, G or V matrix in a mixed or GLS model

Description

These functions call getVarCov in the nlme package. They are intended to have names and functions that are easy to remember.

Usage

getG(fit, ...)

## Default S3 method:
getG(fit, ...)

## S3 method for class 'lme'
getG(fit, ...)

## S3 method for class 'gls'
getG(fit, ...)

getR(fit, ...)

## Default S3 method:
getR(fit, ...)

## S3 method for class 'lme'
getR(fit, ...)

## S3 method for class 'gls'
getR(fit, ...)

getV(fit, ...)

## Default S3 method:
getV(fit, ...)

## S3 method for class 'lme'
getV(fit, ...)

## S3 method for class 'gls'
getV(fit, ...)

Arguments

fit	model created with lme or gls
individuals	for getR or getV with lme objects, select the clusters for which variances are returned. If not specified, the variance of the first cluster is returned.

Value

For lme objects, getG returns the between-cluster variance of random effects, getV, and getR returns a list with the within-cluster marginal variance and the within-cluster conditional variance respectively for the the clusters listed in individuals. If individuals is missing, the variance of the first cluster is returned. ISSUE: For gls objects all functions return the same thing but uninformatively if correlation is clustered and if weights produce differenct variances in the corresponding positions in different clusters.

Methods (by class)

• default: default method

• lme: lme method

• gls: gls method

• default: default method

• lme: lme method

• gls: gls method

• default: gls method

• lme: lme method

• gls: gls method

Examples

library(spida2)
library(nlme)
library(gnew)
data <- expand.grid( Xdev = c(-3,-2,-1,0,1,2,3), id = 1:5 )
 
set.seed(12345)
data <- within(data, {
  Xmean <- 2*id
  X <- Xdev + Xmean
  Y <- (-1 + .1*rnorm(max(id)))[id] * Xdev + 
    2 * Xmean + .3 * id * rnorm(length(id))
})

library(lattice)
gd()
xyplot(Y ~ X, data, groups = id)
fit0 <- lme(Y ~ X, data,
            random = ~ 1+ X |id)
fit <- lme(Y ~ X, data, 
           random = ~ 1 + X | id,
           weights = varConstPower(form = ~ fitted(.)),
           correlation = corAR1(form = ~ 1 | id),
           control = list(returnObject = TRUE))
fitgls <- gls(Y ~ X, data, weights = varConstPower(form = ~ fitted(.)),
           correlation = corAR1(form = ~ 1|id),
           control = list(returnObject = TRUE, maxIter= 1000, 
           verbose = TRUE, msMaxIter = 1000,
               msVerbose=TRUE))
summary(fit)
getVarCov(fit)
getVarCov(fit, individuals = '2')
getVarCov(fit, individuals = '2', type = 'conditional') %>% 
  .[[1]] %>% 
  diag
getVarCov(fit,  type = 'conditional')%>% 
  .[[1]] %>% 
  diag

getG(fit)
getR(fit)[[1]]
getV(fit)[[1]]


(Z <- cbind(1, 2+seq(-3,3)))
Z
(getG(fit))

Z %*% getG(fit) %*% t(Z)

getV(fit)[[1]]
getR(fit)[[1]]
sigma(fit)
getVarCov(fit, type = 'random.effects')
getVarCov(fit)
Z %*% getG(fit) %*% t(Z)
getV(fit)[[1]] - Z %*% getG(fit) %*% t(Z) - getR(fit)[[1]]

getG(fit0)
Z %*% getG(fit0) %*% t(Z)
Z %*% getG(fit0) %*% t(Z) %>% svd %>% .$d
getR(fit0)
sigma(fit0)
getV(fit0)
Z %*% getG(fit0) %*% t(Z) + getR(fit0)[[1]]

gmonette/gnew documentation built on July 9, 2022, 12:57 p.m.