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

In RMark: R Code for Mark Analysis

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
)

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

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.

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

Author(s)

Jeff Laake

Examples

# 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)

RMark documentation built on Aug. 14, 2022, 1:05 a.m.