RMark (version 3.0.0)

var.components.reml: Variance components estimation using REML or maximum likelihood

Description

Computes estimated effects, standard errors and variance components for a set of estimates

Usage

var.components.reml(
  theta,
  design,
  vcv = NULL,
  rdesign = NULL,
  initial = NULL,
  interval = c(-25, 10),
  REML = TRUE
)

Value

A list with the following elements

neglnl

negative log-likelihood for fitted model

AICc

small sample corrected AIC for model selection

sigma

variance component estimates; if rdesign=NULL, only an iid error; otherwise, iid error and random effect error

beta

dataframe with estimates and standard errors of betas for design

vcv.beta

variance-covariance matrix for beta

Arguments

theta

vector of parameter estimates

design

design matrix for fixed effects combining parameter estimates

vcv

estimated variance-covariance matrix for parameters

rdesign

design matrix for random effect (do not use intercept form; eg use ~-1+year instead of ~year); if NULL fits only iid error

initial

initial values for variance components

interval

interval bounds for log(sigma) to help optimization from going awry

REML

if TRUE uses reml else maximum likelihood

Author

Jeff Laake

Details

The function var.components uses method of moments to estimate a single process variance but cannot fit a more complex example. It can only estimate an iid process variance. However, if you have a more complicated structure in which you have random year effects and want to estimate a fixed age effect then var.components will not work because it will assume an iid error rather than allowing a common error for each year as well as an iid error. This function uses restricted maximum likelihood (reml) or maximum likelihood to fit a fixed effects model with an optional random effects structure. The example below provides an illustration as to how this can be useful.

Examples

# \donttest{
# This example is excluded from testing to reduce package check time
# Use dipper data with an age (0,1+)/time model for Phi
data(dipper)
dipper.proc=process.data(dipper,model="CJS")
dipper.ddl=make.design.data(dipper.proc,
   parameters=list(Phi=list(age.bins=c(0,.5,6))))
levels(dipper.ddl$Phi$age)=c("age0","age1+")
md=mark(dipper,model.parameters=list(Phi=list(formula=~time+age)),delete=TRUE)
# extract the estimates of Phi 
zz=get.real(md,"Phi",vcv=TRUE)
# assign age to use same intervals as these are not copied 
# across into the dataframe from get.real
zz$estimates$age=cut(zz$estimates$Age,c(0,.5,6),include=TRUE)
levels(zz$estimates$age)=c("age0","age1+")
z=zz$estimates
# Fit age fixed effects with random year component and an iid error
var.components.reml(z$estimate,design=model.matrix(~-1+age,z),
        zz$vcv,rdesign=model.matrix(~-1+time,z)) 
# Fitted model assuming no covariance structure to compare to 
# results with lme
xx=var.components.reml(z$estimate,design=model.matrix(~-1+age,z),
 matrix(0,nrow=nrow(zz$vcv),ncol=ncol(zz$vcv)),
 rdesign=model.matrix(~-1+time,z)) 
xx
sqrt(xx$sigmasq)
library(nlme)
nlme::lme(estimate~-1+age,data=z,random=~1|time)
# }