Nothing
#' Mixed effect model contstruction
#'
#' Functions that develop structures needed for a mixed effect model
#'
#' mixed.model.admb - creates design matrices and supporting index matrices
#' for use of mixed model in ADMB
#'
#' mixed.model - creates design matrices and supporting index matrices
#' in an alternate list format that is not as easily used in ADMB
#'
#' mixed.model.dat - writes to data file (con) for fixed and random effect stuctures
#'
#' reindex - creates indices for random effects that are specific to the individual capture
#' history; it takes re.indices, splits them by id and creates
#' a ragged array by id (used.indices) with the unique values for that id. index.counts is the number
#' of indices per id to read in ragged array. It then changes re.indices to be an index
#' to the indices within the id from 1 to the number of indices within the id.
#'
#' @usage mixed.model.admb(formula,data)
#'
#' mixed.model(formula,data,indices=FALSE)
#'
#' mixed.model.dat(x,con,idonly,n)
#'
#' reindex(x,id)
#'
#' @aliases mixed.model.admb mixed.model mixed.model.dat reindex
#' @param formula formula for mixed effect mode in the form used in lme4; ~fixed +(re1|g1) +...+(ren|gn)
#' @param data dataframe used to construct the design matrices from the formula
#' @param x list structure created by mixed.model.admb
#' @param con connection to data file which contents will be appended
#' @param id vector of factor values used to split the data up by individual capture history
#' @param idonly TRUE, if random effects not crossed
#' @param n number of capture history records
#' @param indices if TRUE, outputs structure with indices into dm for random effects
#' @return mixed.model.admb returns a list with elements re.dm, a combined design matrix for all of the random effects; and
#' re.indices, matrix of indices into a single vector of random effects to be applied to the
#' design matrix location.
#' mixed.model returns a list (re.list) with an element for each random effect structure. The contents
#' are a standard design matrix (re.dm) if indices==FALSE and a re.dm and re.indices which matches the
#' structure of mixed.model.admb. mixed.model will be more useful with R than ADMB.
#' @author Jeff Laake
#'
mixed.model.admb=function(formula,data)
{
# parse formula for fixed and random effects
mlist=proc.form(formula)
# construct design matrix for fixed effects
# fixed.dm=model.matrix(as.formula(mlist$fix.model),data)
# remainder of code is for random effects unless NULL
reindex=0
re.dm=NULL
re.indices=NULL
freq=NULL
sigma_names=NULL
if(!is.null(mlist$re.model))
{
re_names = sub("^\\s+", "", sapply(strsplit(names(mlist$re.model),"\\|"), function(x) x[2]))
# Loop over each random effect component
for(i in 1:length(mlist$re.model))
{
# Make sure each variable used to define random effect group is a factor variable
if(!all(sapply(model.frame(as.formula(mlist$re.model[[i]]$sub),data),is.factor)))
warning(paste("\n one or more variables in",mlist$re.model[[i]]$sub,"is not a factor variable\n"))
# else
# {
# Compute design matrix for grouping variables
zz=model.matrix(as.formula(mlist$re.model[[i]]$sub),data)
# Not all combinations of factor variable(s) may be used so only use those observed in the data
used.columns=which(colSums(zz)>0)
nre=length(used.columns)
# Compute the indices for this particular grouping structure and reindex if any missing
indices=rowSums(zz*col(zz))
if(nre!=ncol(zz))indices=match(indices,used.columns)
# Compute the design matrix for the random effect formula
zz=model.matrix(as.formula(mlist$re.model[[i]]$model),data)
# Now shift indices to refer to a single vector of random effects across all re groupings
ng=max(indices)
indices=matrix(indices,nrow=length(indices),ncol=ncol(zz))+reindex
indices=t(t(indices)+cumsum(c(0,rep(ng,ncol(zz)-1))))
reindex=max(indices)
# Bind random effect design matrices (re.dm), indices into random effects vector (re.indices) and
# index for the random effect sigma parameter (re.sigma)
colnames(zz)=paste(re_names[i],colnames(zz),sep=":")
re.dm=cbind(re.dm,zz)
re.indices=cbind(re.indices,indices)
# }
}
crossed=sapply(1:ncol(re.indices),function(i) all(table(data$id,re.indices[,i])==1))
ff=re.dm
ff[ff!=0]=1
ff[,crossed]=0
ff=ff*data$freq
zz=unique(cbind(as.vector(re.indices),as.vector(ff)))
zz[zz==0]=1
freq=zz[,2]
}
return(list(re.dm=re.dm,re.indices=re.indices,freq=freq))
# return(list(fixed.dm=fixed.dm,re.dm=re.dm,re.indices=re.indices))
}
mixed.model=function(formula,data,indices=FALSE)
{
mlist=proc.form(formula)
# fixed.dm=model.matrix(as.formula(mlist$fix.model),data)
re.list=NULL
if(!is.null(mlist$re.model))
{
re.list=vector("list",length=length(mlist$re.model))
for(i in 1:length(mlist$re.model))
{
if(!all(sapply(model.frame(as.formula(mlist$re.model[[i]]$sub),data),is.factor)))
warning(paste("\n one or more variables in",mlist$re.model[[i]]$sub,"is not a factor variable\n"))
# else
# {
if(indices)
{
zz=model.matrix(as.formula(mlist$re.model[[i]]$sub),data)
used.columns=which(colSums(zz)>0)
nre=length(used.columns)
indices=rowSums(zz*col(zz))
if(nre!=ncol(zz))indices=match(indices,used.columns)
re.dm=model.matrix(as.formula(mlist$re.model[[i]]$model),data)
re.list[[i]]=list(re.indices=indices,re.dm=re.dm)
names(re.list)=names(mlist$re.model)
}else
{
sub=strsplit(mlist$re.model[[i]]$sub, "~ ")[[1]][2]
mod=strsplit(mlist$re.model[[i]]$model, "~ ")[[1]][2]
# strip leading and trailing blank space
mod=gsub("^\\s+|\\s+$", "", mod)
if(mod=="1") form=formula(paste(c("~", sub), collapse=""))
else form=formula(paste(c("~", mod, ":(", sub, ")"), collapse=""))
dm=model.matrix(form, data)
colnames(dm) = paste("(",mod,"|",strsplit(sub, " - 1"),")",c(1:ncol(dm)), sep="")
dm=dm[,colSums(dm)!=0]
re.list[[i]]=Matrix(dm)
names(re.list)[i]=paste(c(mod, "|",strsplit(sub, " - 1")), collapse="")
}
# }
}
}
return(list(re.list=re.list))
}
mixed.model.dat=function(x,con,idonly,n)
{
# number of columns of fixed dm
#write(ncol(x$fixed.dm),con,append=TRUE)
# fixed dm
#write(t(x$fixed.dm),con,ncolumns=ncol(x$fixed.dm),append=TRUE)
if(!is.null(x$re.dm))
{
# number of random effects
if(is.null(x$index.counts))
if(!idonly)
write(max(x$re.indices),con,append=TRUE)
else
write(n,con,append=TRUE)
else
write(max(sapply(x$used.indices,function(x) ifelse(length(x)>0,max(x),-Inf))),con,append=TRUE)
# number of columns of re dm
write(ncol(x$re.dm),con,append=TRUE)
# re dm
write(t(x$re.dm),con,ncolumns=ncol(x$re.dm),append=TRUE)
if(!idonly)
{
# re indices
write(t(x$re.indices),con,ncolumns=ncol(x$re.indices),append=TRUE)
# index counts and used.indices, if not null
if(!is.null(x$index.counts))
{
write(x$index.counts,con,append=TRUE)
for(i in 1:length(x$used.indices))
if(length(x$used.indices[[i]])>0)
write(x$used.indices[[i]],con,ncolumns=length(x$used.indices[[i]]),append=TRUE)
}
}
}
else
{
# number of re =0
write(0,con,append=TRUE)
# number of columns of re dm=0
write(0,con,append=TRUE)
}
invisible()
}
reindex=function(x,id)
{
# accept indices; return list with new indices and
vec.indices=as.vector(x$re.indices)
# re-index
xx=t(apply(x$re.indices,1,match,table=sort(unique(vec.indices[!is.na(vec.indices)]))))
if(nrow(xx)==1)xx=t(xx)
# loop over individuals getting unique non-NA indices; save those in ragged array
x$used.indices=lapply(split(xx,id),function(z) sort(unique(z[!is.na(z)])))
x$index.counts=sapply(x$used.indices,length)
# re-index values into unique indices
sp=split(1:nrow(xx),id)
final=vector("list",length(sp))
for (i in 1:length(sp))
final[[i]]=t(apply(xx[sp[[i]],,drop=FALSE],1,function(z) match(z,x$used.indices[[i]])))
x$re.indices=do.call("rbind",final)
x$re.indices[is.na(x$re.indices)]=0
return(x)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.