var.components.reml {RMark} | R Documentation |
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)