Package 'RMark'

Description

The RMark package is a collection of R functions that can be used as an interface to MARK for analysis of capture-recapture data. Some initial ideas are described here but you are strongly encouraged to read use the help files with this package and to read the documentation at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip The above link appears the first time you enter library(RMark) in an R session.

Details

The RMark package contains various functions that import/export capture data, build capture-recapture models, run the FORTRAN program MARK.EXE, and extract and display output. Program MARK has its own user interface; however, model development can be rather tedious and error-prone because the parameter structure and design matrix are created by hand. This interface in R was created to use the formula and design matrix functions in R to ease model development and reduce errors. This R interface has the following advantages: 1) Uses model notation to create design matrices rather than designing them by hand in MARK or in EXCEL, which makes model development faster and more reliable. All-different PIMS are automatically created for each group (if any). 2) Allows models based on group (factor variables) and individual covariates with groups created on the fly. Age, cohort, group and time variables are pre-defined for use in formulas. 3) Both real and beta labels are automatically added for easy output interpretation. 4) Input, output and specific results (eg parameter estimates, AICc etc) are stored in an R object where they can be manipulated as deemed useful (eg plotting, further calculations, simulation etc). 5) Parameter estimates can be displayed in triangular PIM format (if appropriate) for ease of interpretation. 6) Easy setup of batch jobs and the calls to the R functions document the model specifications and allow models to be easily reproduced or re-run if data are changed. 7) Covariate-specific estimates of real parameters can be computed within R without re-running the analysis.

The MARK capture-recapture models that are currently supported are provided in MarkModels.pdf which is installed in the RMark directory of your R library. You can also find a list in MARK under Help/Data Types. ' There is one limitation of this interface. All models in this interface are developed via a design matrix approach rather than coding the model structure via parameter index matrices (PIMS). In most cases, a logit or other link is used by default which has implications for ability of MARK to count the number of identifiable parameters (see dipper for an example). However, beginning with v1.7.6 the sin link is now supported if the formula specifies an identity design matrix for the parameter.

Before you begin, you must have installed MARK (http://www.phidot.org/software/mark/) on your computer or at least have a current copy of MARK.EXE. As long as you selected the default location for the MARK install (c:/Program Files/Mark), the RMark library will be able to find it. If for some reason, you choose to install it in a different location, see the note section in mark for instructions on setting the variable MarkPath to specify the path. In addition to installing MARK, you must have installed the RMark library into the R library directory. Once done with those tasks, run R and enter library(RMark) (or put it in your .First function) to attach the library of functions.

The following is a categorical listing of the functions in the package with a link to the help for each function. To start, read the help for functions import.chdata and mark to learn how to import your data and fit a simple model. The text files for the examples shown in import.chdata are in the subdirectory data within the R Library directory in RMark. Next look at the example data sets and analyses dipper, edwards.eberhardt, and example.data. After you see the structure of the examples and the use of functions to fit a series of analyses, explore the remaining functions under Model Fitting, Batch Analyses, Model Selection and Summary and Display. If your data and models contain individual covariates, read the section on Real Parameter Computation to learn how to compute estimates of real parameters at various covariate values.

Input/Output data & results

import.chdata,read.mark.binary, extract.mark.output

Exporting Models to MARK interface

export.MARK

Model Fitting

mark, process.data, make.design.data, add.design.data, make.mark.model, run.mark.model merge_design.covariates

Batch analyses with functions

run.models, collect.models, create.model.list, mark.wrapper

Summary and display

summary.mark, print.mark, print.marklist, get.real, compute.real, print.summary.mark

Model Selection/Goodness of fit

adjust.chat, adjust.parameter.count, model.table , release.gof, model.average

Real Parameter computation

find.covariates, fill.covariates, compute.real , covariate.predictions

Utility and internal functions

mark.wrapper,collect.model.names, compute.design.data, extract.mark.output, inverse.link, deriv.inverse.link, setup.model, setup.parameters, valid.parameters, cleanup,TransitionMatrix,

For examples, see dipper for CJS and POPAN, see example.data for CJS with multiple grouping variables, see edwards.eberhardt for various closed-capture models, see mstrata for Multistrata, and see Blackduck for known fate. The latter two are examples of the use of mark.wrapper for a shortcut approach to creating a series of models. Other examples have been added for the various other models. In MarkModels.pdf it also lists the name of examples that are provided for each model.

Author(s)

Jeff Laake

References

MARK: Dr. Gary White, Department of Fishery and Wildlife Biology, Colorado State University, Fort Collins, Colorado, USA http://www.phidot.org/software/mark/

Description

Creates new design data fields in (ddl) that bin the fields cohort, age or time. Other fields (e.g., effort value for time) can be added to ddl with R commands. Also see merge_design.covariates to merge data like environmental variables or any other data that is common to all animals or group/time specific values.

Usage

add.design.data(
  data,
  ddl,
  parameter,
  type = "age",
  bins = NULL,
  name = NULL,
  replace = FALSE,
  right = TRUE
)

Arguments

Details

Design data can be added to the parameter specific design dataframes with R commands and this function does NOT have to be used to add design data. However, often the additional fields will be functions of cohort, age or time. add.design.data provides an easy way to add fields that bin (put into intervals) the original values of cohort, age or time. For example, age may have levels from 0 to 10 which means the formula ~age will have 11 parameters, one for each level of the factor. It might be more desirable and more parsimonious to have a simpler 2 age class model of young and adults. This can be done easily by adding a new design data field that bins age into 2 intervals (age 0 and 1+) as in the following example:

 ddl=make.design.data(proc.example.data)
ddl=add.design.data(proc.example.data,ddl,parameter="Phi",type="age",
bins=c(0,.5,10),name="2ages") 

By default, the bins are open on the left and closed on the right (i.e., binning x by (x1,x2] is equivalent to x1<x<=x2) except for the first interval which is closed on the left. Thus, for the above example, the age bins are [0,.5] and (.5,10]. Since the ages in the example are 0,1,2... using any value >0 and <1 in place of 0.5 would bin the ages into 2 classes of 0 and 1+. This behavior can be modified by changing the argument right=FALSE to create an interval that is closed on the left and open on the right. In some cases this can make reading the values of the levels somewhat easier. It is important to recognize that the new variable is only added to the design data for the defined parameter and can only be used in model formula for that parameter. Multiple calls to add.design.data can be used to add the same or different fields for the various parameters in the model. For example, the same 2 age class variable can be added to the design data for p with the command:

ddl=add.design.data(proc.example.data,ddl,parameter="p",type="age",
bins=c(0,.5,10),name="2ages") 

The name must be unique within the parameter design data, so they should not use pre-defined values of group, age, Age, time, Time, cohort, Cohort. If you choose a name that already exists in the design data for the parameter, it will not be added but it can replace the variable if replace=TRUE. For example, the 2ages variable can be re-defined to use 0-1 and 2+ with the command:

ddl=add.design.data(proc.example.data,ddl,parameter="Phi",type="age",
bins=c(0,1,10),name="2ages",replace=TRUE) 

Keep in mind that design data are stored with the mark model object so if a variable is redefined, as above, this could become confusing if some models have already been constructed using a different definition for the variable. The model formula and names would appear to be identical but they would have a different model structure. The difference would be apparent if you examined the design data and design matrix of the model object but would the difference would be transparent based on the model names and formula. Thus, it would be best to avoid constructing models from design data fields with different structures but the same name.

Value

Design data list with new field added for the specified parameter. See make.design.data for a description of the list structure.

Note

For the specific case of "closed" capture models, the parameters p (capture probability) and c (recapture probability) can be treated in a special fashion. Because they really the same type of parameter, it is useful to be able to share a common model structure (i.e., same columns in the design matrix). This is indicated with the share=TRUE element in the model description for p. If the parameters are shared then the additional covariate c is added to the design data, which is c=0 for parameter p and c=1 for parameter c. This enables an additive model to be developed where recapture probabilities mimic the pattern in capture probabilities except for an additive constant. The covariate c can only be used in the model for p if share=TRUE. If the latter is not set using c in a formula will result in an error. Likewise, if share=TRUE, then the design data for p and c must be the same because the design data are merged in constructing the design matrix. Thus if you add design data for parameter p, you should add a similar field for parameter c if you intend to fit shared models for the two parameters. If the design data do not match and you try to fit a shared model, an error will result.

Author(s)

Jeff Laake

See Also

make.design.data, process.data

Examples

# This example is excluded from testing to reduce package check time
data(example.data)
example.data.proc=process.data(example.data)
ddl=make.design.data(example.data.proc)
ddl=add.design.data(example.data.proc,ddl,parameter="Phi",type="age",
  bins=c(0,.5,10),name="2ages")
ddl=add.design.data(example.data.proc,ddl,parameter="p",type="age",
bins=c(0,.5,10),name="2ages")
ddl=add.design.data(example.data.proc,ddl,parameter="Phi",type="age",
bins=c(0,1,10),name="2ages",replace=TRUE)

Description

Modifies number of estimated parameters and the resulting AICc value for model selection.

Usage

adjust.parameter.count(model, npar)

Arguments

Details

When a model is run the parameter count determined by MARK is stored in results$npar and the AICc value is stored in results$AICc. If the argument adjust is set to TRUE in the call to run.mark.model and MARK determined that the design matrix was not full rank (i.e., the parameter count is less than the columns of the design matrix), then the parameter count from MARK is stored in results$npar.unadjusted and AICc in results$AICc.unadjusted and results$npar is set to the number of columns of the design matrix and results$AICc uses the assumed full rank value of npar. This function allows the parameter count to be reset to any value less than or equal to the number of columns in the design matrix. If results$npar.unadjusted exists it is kept as is. If it doesn't exist, then the current values of results$npar and results$AICc are stored in the .unadjusted fields to maintain the values from MARK, and the new adjusted values defined by the function argument npar are stored in results$npar and results$AICc. In the example below, the CJS model Phi(t)p(t) is fitted with the call to mark which defaults to adjust=TRUE. This is used to show how adjust.parameter.count can be used to adjust the count to 11 from the full rank count of 12. Alternatively, the argument adjust=FALSE can be added to prevent the adjustment which is appropriate in this case because Phi(6) and p(6) are confounded.

Value

model: the mark model object with the adjustments made

Author(s)

Jeff Laake

See Also

run.mark.model,model.table

Examples

# This example is excluded from testing to reduce package check time
data(dipper)
ptime=list(formula=~time)
Phitime=list(formula=~time)
dipper.phitime.ptime=mark(dipper,model.parameters=list(Phi=Phitime, p=ptime),delete=TRUE)
dipper.phitime.ptime=adjust.parameter.count(dipper.phitime.ptime,11)
dipper.phitime.ptime=mark(dipper,model.parameters=list(Phi=Phitime, p=ptime),
                           adjust=FALSE,delete=TRUE)

Description

Adjust value of over-dispersion constant or another result value for a collection of models which modifies model selection criterion and estimated standard errors.

Usage

adjust.value(field="n",value,model.list)
       adjust.chat(chat=1,model.list)

Arguments

Details

The value of chat is stored with the model object except when there is no over-dispersion (chat=1). This function assigns a new value of chat for the collection of models specified by model.list and/or type. The value of chat is used by model.table for model selection in computing QAICc unless chat=1. It is also used in summary.mark, get.real and compute.real to adjust standard errors and confidence intervals. Note that the standard errors and confidence intervals in results$beta,results$beta.vcv results$real, results$derived and results$derived.vcv are not modified and always assume chat=1.

It can also be used to modify a field in model$results such as n which is ESS (effective sample size) from MARK output that is used in AICc/QAICc calculations.

Value

model.list with all models given the new chat value and model.table adjusted for chat values

Note

See note in collect.models

Author(s)

Jeff Laake

See Also

model.table, summary.mark, get.real ,compute.real

Examples

#
# The following are examples only to demonstrate selecting different 
# model sets for adjusting chat and showing model selection table. 
# It is not a realistic analysis.
#

# This example is excluded from testing to reduce package check time
data(dipper)
do_example=function()
{
mod1=mark(dipper,delete=TRUE)
mod2=mark(dipper,model.parameters=list(Phi=list(formula=~time)),delete=TRUE)
mod3=mark(dipper,model="POPAN",initial=1,delete=TRUE)
cjs.results=collect.models(type="CJS")
cjs.results  # show model selection results for "CJS" models
}
cjs.results=do_example()
cjs.results
# adjust chat for all models to 2
cjs.results=adjust.chat(2,cjs.results) 
cjs.results

Description

A known fate data set on Black ducks that accompanies MARK as an example analysis using the Known model.

Format

A data frame with 48 observations on the following 5 variables.

ch

a character vector containing the encounter history of each bird

BirdAge

the age of the bird: a factor with levels 0 1 for young and adult

Weight

the weight of the bird at time of marking

Wing_Len

the wing-length of the bird at time of marking

condix

the condition index of the bird at time of marking

Details

This is a data set that accompanies program MARK as an example for Known fate. The data can be stratified using BirdAge as a grouping variable. The function run.Blackduck defined below in the examples creates some of the models used in the dbf file that accompanies MARK.

Note that in the MARK example the variable is named Age. In the R code, the fields "age" and "Age" have specific meanings in the design data related to time since release. These will override the use of a field with the same name in the individual covariate data, so the names "time", "Time", "cohort", "Cohort", "age", and "Age" should not be used in the individual covariate data with possibly the exception of "cohort" which is not defined for models with "Square" PIMS such as POPAN and other Jolly-Seber type models.

Examples

data(Blackduck)
# Change BirdAge to numeric; starting with version 1.6.3 factor variables are
# no longer allowed.  They can work as in this example but they can be misleading
# and fail if the levels are non-numeric.  The real parameters will remain 
# unchanged but the betas will be different.
Blackduck$BirdAge=as.numeric(Blackduck$BirdAge)-1
run.Blackduck=function()
{
#
# Process data
#
bduck.processed=process.data(Blackduck,model="Known")
#
# Create default design data
#
bduck.ddl=make.design.data(bduck.processed)
#
#  Add occasion specific data min < 0; I have no idea what it is
#
bduck.ddl$S$min=c(4,6,7,7,7,6,5,5)
#
#  Define range of models for S
#
S.dot=list(formula=~1)
S.time=list(formula=~time)
S.min=list(formula=~min)
S.BirdAge=list(formula=~BirdAge)
#
# Note that in the following model in the MARK example, the covariates
# have been standardized.  That means that the beta parameters will be different
# for BirdAge, Weight and their interaction but the likelihood and real parameter
# estimates are the same.
#
S.BirdAgexWeight.min=list(formula=~min+BirdAge*Weight)
S.BirdAge.Weight=list(formula=~BirdAge+Weight)
#
# Create model list and run assortment of models
#
model.list=create.model.list("Known")
bduck.results=mark.wrapper(model.list,data=bduck.processed,ddl=bduck.ddl,
               invisible=FALSE,threads=1,delete=TRUE)

#
# Return model table and list of models
#
return(bduck.results)
}
bduck.results=run.Blackduck()
bduck.results

Description

A recovery data set for mallards in San Luis Valley Colorado

Format

A data frame with 108 observations on the following 5 variables.

ch

a character vector containing the encounter history of each bird

freq

frequency of that encounter history

ReleaseAge

the age of the bird when it was released

Details

This is a data set that accompanies program MARK as an example for Brownie and Seber recovery model. In those input files it is in a summarized format. Here it is in the LD encounter history format. The data can be stratified using ReleaseAge (Adult, Young) as a grouping variable.

Note that in the MARK example the variable is named Age. In the R code, the fields "age" and "Age" have specific meanings in the design data related to time since release. These will override the use of a field with the same name in the individual covariate data, so the names "time", "Time", "cohort", "Cohort", "age", and "Age" should not be used in the individual covariate data with possibly the exception of "cohort" which is not defined for models with "Square" PIMS such as POPAN and other Jolly-Seber type models.

Examples

# brownie=import.chdata("brownie.inp",field.types=c("n","f"))

# This example is excluded from testing to reduce package check time
 data(brownie)
# default ordering of ReleaseAge is alphabetic so it is 
# Adult, Young which is why initial.ages=c(1,0)
# Seber Recovery
br=process.data(brownie,model="Recovery",groups="ReleaseAge",age.var=1,initial.ages=c(1,0))
br.ddl=make.design.data(br,parameters=list(S=list(age.bins=c(0,1,10)),
                                           r=list(age.bins=c(0,1,10))),right=FALSE)
mod=mark(br,br.ddl,model.parameters=list(S=list(formula=~-1+age:time,link="sin"),
                                           r=list(formula=~-1+age:time,link="sin")),delete=TRUE)
# Brownie Recovery
br=process.data(brownie,model="Brownie",groups="ReleaseAge",age.var=1,initial.ages=c(1,0))
br.ddl=make.design.data(br,parameters=list(S=list(age.bins=c(0,1,10)),
                               f=list(age.bins=c(0,1,10))),right=FALSE)
mod=mark(br,br.ddl,model.parameters=list(S=list(formula=~-1+age:time,link="sin"),
                               f=list(formula=~-1+age:time,link="sin")),delete=TRUE)
mod=mark(br,br.ddl,model.parameters=list(S=list(formula=~-1+age,link="sin"),
                               f=list(formula=~-1+age,link="sin")),delete=TRUE)
#Random effects Seber recovery (not run to save computation time)
#br=process.data(brownie,model="REDead",groups="ReleaseAge",age.var=1,initial.ages=c(1,0))
#br.ddl=make.design.data(br,parameters=list(S=list(age.bins=c(0,1,10)),
#                                       r=list(age.bins=c(0,1,10))),right=FALSE)
#mod=mark(br,br.ddl,model.parameters=list(S=list(formula=~age),r=list(formula=~age)),delete=TRUE)
#Pledger Mixture Seber recovery
br=process.data(brownie,model="PMDead",groups="ReleaseAge",
                           mixtures=3,age.var=1,initial.ages=c(1,0))
br.ddl=make.design.data(br,parameters=list(S=list(age.bins=c(0,1,10)),
                            r=list(age.bins=c(0,1,10))),right=FALSE)
mod=mark(br,br.ddl,model.parameters=list(pi=list(formula=~mixture),
                     S=list(formula=~age+mixture),r=list(formula=~age)),delete=TRUE)
br=process.data(brownie,model="PMDead",groups="ReleaseAge",
                     mixtures=2,age.var=1,initial.ages=c(1,0))
br.ddl=make.design.data(br,parameters=list(S=list(age.bins=c(0,1,10)),
                      r=list(age.bins=c(0,1,10))),right=FALSE)
mod=mark(br,br.ddl,model.parameters=list(pi=list(formula=~age),
                      S=list(formula=~age+mixture),r=list(formula=~age)),delete=TRUE)

Description

An example of the Burnham live-dead model using simulated data LD1.inp from Chapter 9 of Cooch and White

Author(s)

Luke Eberhart-Phillips<luke.eberhart at gmail.com>

Examples

###############################################################################
#### RMARK script for conducting the Burnham model tutorial in Chapter 9.3 ####
#################### the of the Cooch and White MARK book #####################
###############################################################################
##################### Code by: Luke Eberhart-Phillips #########################
###### Dept. Animal Behaviour, Bielefeld University, Bielefeld, Germany #######
##################### email: luke.eberhart at gmail.com ##########################
###############################################################################

# import/convert the simulated "LD1.inp" MARK capture history into an RMARK 
# dataframe, while defining the two groups as "Y" for individuals marked as 
# young, and "A" for individuals marked as adults
# NOTE: the "LD1.inp" file is found in the zipped folder downloaded when you
# click on "Example data files" in the drop-down menu of the MARK book webpage
#  \url{http://www.phidot.org/software/mark/docs/book/}
pathtodata=paste(path.package("RMark"),"extdata",sep="/")
LD=convert.inp(paste(pathtodata,"ld1",sep="/"),
           group.df=data.frame(age_marked=c("Y","A")))
# process the data by defining the model type as "Burnham" and the groups in
# the data.  In this case the only group is the age at which individuals were 
# marked
LD.proc=process.data(data = LD, 
	model = "Burnham",
	groups=c("age_marked"),
	age.var=1,
	initial.age=c(1,0))

# make the design data from the process data above
LD.ddl=make.design.data(LD.proc)

# add the correct binning to the design data so that individuals that were
# marked as young are adults in their second year of life, where as those
# marked as adults are adults for their entire life.
LD.ddl=add.design.data(data = LD.proc,
	ddl = LD.ddl,
	parameter="S", 
	type = "age",
	bins = c(0,1,8),
	right = FALSE,
	name = "age",
	replace = TRUE)

# do the same to the F parameter
LD.ddl=add.design.data(data = LD.proc,
	ddl = LD.ddl,
	parameter="F", 
	type = "age",
	bins = c(0,1,8),
	right = FALSE,
	name = "age",
	replace = TRUE)

# check parameter matrix to see if groups were binned correctly in the S matrix
PIMS(mark(data = LD.proc,
			ddl = LD.ddl,
			model.parameters=list(S=list(formula=~age)),
			delete=TRUE,
			model = "Burnham"),
	"S")

# Create the formulas that describe variation in the parameter we want to test.
# In this case we want to test for an age effect on survival and fidelity,
# while keeping recapture and recovery probabilities constant.
S.age=list(formula=~age) # S(age)
p.dot=list(formula=~1) # p(.)
F.age=list(formula=~age) # F(age)
r.dot=list(formula=~1) # r(.)

# Run the model
LD.model.age.F.S=mark(data = LD.proc,
	ddl = LD.ddl, 
	model.parameters = list(S = S.age, p = p.dot, 
			F =F.age, r = r.dot), 
	invisible = FALSE, 
	model = "Burnham",delete=TRUE)

# Check the parameter estimates, they should be the same as those generated
# when doing the tutorial in chapter 9.3 of the in MARK Book (table on pg 9-8)
LD.model.age.F.S$results$real

# Clean your working directory
cleanup(ask=FALSE)

Description

Identifies all unused (orphaned) mark*.inp, .vcv, .res and .out and .tmp files in the working directory and removes them. The orphaned files are determined by examining all mark objects and lists of mark objects (created by collect.models) to create a list of files in use. All other files are treated as orphans to delete.

Usage

cleanup(lx = NULL, ask = TRUE, prefix = "mark")

Arguments

Details

This function removes orphaned output files from MARK. This occurs when there are output files in the subdirectory that are not associated with a mark object in the current R session (.Rdata). For example, if you repeat an analysis or set of analyses and store them in the same object then the original set of output files would no longer be linked to an R object and would be orphaned.

As an example, consider the mallard analysis. The first time you run the analysis script in an empty subdirectory it would create 9 sets of MARK output files (mark001.out,.vcv,.res,.inp to mark009.out,.vcv,.res,.inp) and each would be linked to one of the objects in nest.results. When the command AgePpnGrass=nest.results$AgePpnGrass was issued, both of those mark objects were linked to the same set of output files. Now if you were to repeat the above commands and re-run the models and stored the results into nest.results again, it would create files with prefixes 10 to 18. Because that would have replaced nest.results, none of the files numbered 1 to 9 would be linked to an object. cleanup(ask=FALSE) automatically removes those orphan files. If you delete all objects in the R session with the command rm(list=ls(all=TRUE)), then subsequently cleanup(ask=FALSE) will delete all MARK output files because all of them will be orphans. Output files can also become orphans if MARK finishes but for some reason R crashes or you forget to save your session before you exit R. Orphan output files can be re-linked to an R object without re-running MARK by re-running the mark function in R and specifing the filename argument to match the base portion of the orphaned output file (eg "mark067"). It will create all of the necessary R objects and then asks if you want to use the existing file. If you respond affirmatively it will link to the orphaned files.

Value

None

Author(s)

Jeff Laake

Description

MARK model beta parameters

Usage

## S3 method for class 'mark'
coef(object,...)

Arguments

Value

A vector or dataframe of beta estimates

Author(s)

Jeff Laake

Description

Either names of all mark model objects (type=NULL) or names of mark model objects of a specific type (type) are extracted from a vector of R objects (lx) that was collected from the parent environment (frame) of the function that calls collect.model.names. Thus, it is two frames back (parent.frame(2)).

Usage

collect.model.names(lx, type = NULL, warning = TRUE)

Arguments

Details

If type=NULL then the names of all objects of class(x)[1]="mark" in lx are returned. If type is specified, then the names of all objects of class(x)=c("mark",type) in lx are returned.

This function was written with the intention that it would be called from other functions ( e.g., collect.models, run.models) but it will work if called directly (e.g., collect.model.names( lx=ls())). While this function returns a vector of model names, collect.models returns a list of model objects. The latter can be used to easily create a list of models created in a function to be used as a return value without listing all the names of the functions. It uses collect.model.names to perform that function.

Value

model.list: a vector of mark model names

Author(s)

Jeff Laake

See Also

collect.models, run.models, model.table

Description

Collects mark models contained in lx of specified type (if any) and returns models in a list with a table of model results if table=TRUE.

Usage

collect.models(
  lx = NULL,
  type = NULL,
  table = TRUE,
  adjust = TRUE,
  external = FALSE
)

Arguments

Details

If lx is NULL a vector of object names in the parent frame is constructed for lx. Within lx all mark model objects (i.e., class(x)[1]=="mark") are returned if type is NULL. If type is specified and is valid, then the names of all mark model objects of the specified type (i.e., class(x) = c("mark",type)) in lx are returned. If table=TRUE a table of model selection results is also included in the returned list.

This function was written to be able to easily collect a series of mark models in a list without specifying the names of each model object. This is useful in constructing a return value of a function that runs a series of models for a particular analysis. For an example see dipper.

Value

model.list: a list of mark models and optionally a table of model results.

Note

This function and others that use it or use collect.model.names to collect a series of models or assign a value to a series of models (e.g., adjust.chat) should be used with a degree of caution. It is important to understand the scope of the collection. If the call to this function is made at the R prompt, then it will collect all models (of a particular type if any) within the current .Rdata file. If the call to this function (or one like it) is called from within a function that runs a series of analyses then the collection is limited to the function frame (i.e., only models defined within the function). Thus, it is wise to either use a different .Rdata file for each data set (e.g., one for dipper, another for edwards.eberhardt, etc) or to run everything within functions as illustrated by dipper or edwards.eberhardt. Using a separate .Rdata file is equivalent to having separate .DBF/.FPT files with MARK. It is important to note that functions such as adjust.chat will adjust the value of chat across analyses unless specifically given a list of models to adjust.

Author(s)

Jeff Laake

See Also

merge.mark,remove.mark, collect.model.names,run.models,model.table,dipper

Examples

# see examples in dipper, edwards.eberhardt and example.data

Description

For a specific type of parameter (e.g., Phi, p, r etc), it creates a data frame containing design data for each parameter of that type in the model as structured by an all different PIM (parameter information matrix). The design data are used in constructing the design matrix for MARK with user-specified model formulae as in make.mark.model.

Usage

compute.design.data(
  data,
  begin,
  num,
  type = "Triang",
  mix = FALSE,
  rows = 0,
  pim.type = "all",
  secondary,
  nstrata = 1,
  tostrata = FALSE,
  strata.labels = NULL,
  subtract.stratum = strata.labels,
  common.zero = FALSE,
  sub.stratum = 0,
  limits = NULL,
  events = NULL,
  use.events = NULL,
  mscale = 1,
  subtract.events = NULL
)

Arguments

Details

This function is called by make.design.data to create all of the default design data for a particular type of model and by add.design.data to add binned design data fields for a particular type of parameter. The design data created by this function include group, age, time and cohort as factors variables and continuous (non-factor) versions of all but group. In addition, if groups have been defined for the data, then a data column is added for each factor variable used to define the groups. Also for specific closed capture heterogeneity models (model="HetClosed", "FullHet", "HetHug", "FullHetHug") the data column mixture is added to the design data. The arguments for this function are defined for each model by the function setup.model.

Value

design.data: a data frame containing all of the design data fields for a particular type of parameter

Author(s)

Jeff Laake

See Also

make.design.data, add.design.data

Description

Computes link values (design*beta) for real parameters, and var-cov from design matrix (design) and coefficients (beta)

Usage

compute.link(
  model,
  beta = NULL,
  design = NULL,
  data = NULL,
  parm.indices = NULL,
  vcv = TRUE
)

Arguments

Details

This function is very similar to compute.real except that it provides estimates of link values before they are transformed to real estimates using the inverse-link. It is called by get.link to make calculations but can be called separately. The value is always a dataframe for the estimates and design data and optionally a variance-covariance matrix. See get.real for further details about the arguments.

Value

estimates: If vcv=TRUE, a list is returned with elements vcv.link and the dataframe estimates. If vcv=FALSE, only the estimates dataframe is returned which has the same structure as in get.real.

Author(s)

Jeff Laake

See Also

get.link

Description

Computes link values from reals using 1-1 real to beta(=link) transformation. Also, creates a v-c matrix for the link values if vcv.real is specified.

Usage

compute.links.from.reals(
  x,
  model,
  parm.indices = NULL,
  vcv.real = NULL,
  use.mlogits = TRUE
)

Arguments

Details

It has 2 uses both related to model averaged estimates. Firstly, it is used to transform model averaged estimates so the normal confidence interval can be constructed on the link values and then back-transformed to real space. The second function is to enable parametric bootstrapping in which the error distbution is assumed to be multivariate normal for the link values. From a single model, the link values are easily constructed from the betas and design matrix so this function is not needed. But for model averaging there is no equivalent because the real parameters are averaged over a variety of models with the same real parameter structure but differing design structures. This function allows for link values and their var-cov matrix to be created from the model averaged real estimates.

Value

A list with the estimates (link values) and the links that were used. If vcv.real = TRUE, then the v-c matrix of the links is also returned.

Author(s)

Jeff Laake

See Also

model.average

Description

Computes real estimates and var-cov from design matrix (design) and coefficients (beta) using specified link functions

Usage

compute.real(
  model,
  beta = NULL,
  design = NULL,
  data = NULL,
  se = TRUE,
  vcv = FALSE
)

Arguments

Details

The estimated real parameters can be derived from the estimated beta parameters, a completed design matrix, and the link function specifications. MARK produces estimates of the real parameters, se and confidence intervals but there are at least 2 situations in which it is useful to be able to compute them after running the analysis in MARK: 1) adjusting confidence intervals for estimated over-dispersion, and 2) making estimates for specific values of covariates. The first case is done in get.real with a call to this function. It is done by adjusting the estimated standard error of the beta parameters by multiplying it by the square root of chat to adjust for over-dispersion. A normal 95 +/- 1.96*se) and this is then back-transformed to the real parameters using inverse.link with the appropriate inverse link function for the parameter to construct a 95 There is one exception. For parameters using the mlogit transformation, a logit transformation of each individual real Psi and its se are used to derive the confidence interval. The estimated standard error for the real parameter is also scaled by the square root of the over-dispersion constant chat stored in model$chat. But, the code actually computes the variance-covariance matrix rather than relying on the values from the MARK output because real estimates will depend on any individual covariate values used in the model which is the second reason for this function.

New values of the real parameter estimates can easily be computed by simply changing the values of the covariate values in the design matrix and computing the inverse-link function using the beta parameter estimates. The covariate values to be used can be specified in one of 2 ways. 1) Prior to making a call to this function, use the functions find.covariates to extract the rows of the design matrix with covariate values and either fill in those values aautomatically with the options provided by find.covariates or edit those values to be the ones you want and then use fill.covariates to replace the values into the design matrix and use it as the value for the argument design, or 2) automate this step by specifying a value for the argument data which is used to take averages of the covariate values to fill in the covariate entries of the design matrix. In computing real parameter estimates from individual covariate values it is important to consider the scale of the individual covariates. By default, an analysis with MARK will standardize covariates by subtracting the mean and dividing by the standard deviation of the covariate value. However, in the RMark library all calls to MARK.EXE do not standardize the covariates and request real parameter estimates based on the mean covariate values. This was done because there are many instances in which it is not wise to use the standardization implemented in MARK and it is easy to perform any standardization of the covariates with R commands prior to fitting the models. Also, with pre-standardized covariates there is no confusion in specifying covariate values for computation of real estimates. If the model contains covariates and the argument design is not specified, the design matrix is extracted from model and all individual covariate values are assigned their mean value to be consistent with the default in the MARK analysis.

If a value for beta is given, those values are used in place of the estimates model$results$beta$estimate.

Value

A data frame (real) is returned if vcv=FALSE; otherwise, a list is returned also containing vcv.real:

Author(s)

Jeff Laake

See Also

get.real,fill.covariates,find.covariates,inverse.link,deriv_inverse.link

Examples

# see examples in fill.covariates

Description

Converts encounter history inp files used to create a MARK project into a dataframe for use with RMark. Group structure in frequencies is converted to factor variables that can be used to create groups in RMark. Covariates are copied straight across. Only works with encounter history format for input files and not specialized ones for known-fate or Brownie models.

Usage

convert.inp(
  inp.filename,
  group.df = NULL,
  covariates = NULL,
  use.comments = FALSE
)

Arguments

Details

The encounter history format for MARK is structured as follows: capture (encounter) history, followed by a frequency field for each group, followed by any covariates and then a semi-colon at the end of the line. Comments are allowed within /* and */. The RMark format is a dataframe with a different structure. Each record(row) in the dataframe is for one or more animals within a single group and if there is group structure then the dataframe contains factor variables that can be used to create groups. For example, the following is a little snippet of the same data with 2 groups Males/Females and a covariate weight in the two different formats:

 MARK encounter history file (in make believe test.inp): 1001
1 0 10; 1101 0 2 5; 0101 3 1 6;

RMark dataframe: ch freq sex weight 1001 1 M 10 1101 2 F 5 0101 3 M 6 0101 1
F 6 

To convert from the MARK format to the RMark format it is necessary to define the variables used to define the groups (if any) and to define the covariate field names (if any). For the example above, if test.inp is in the same directory as the current working directory, the call would be:

test = convert.inp("test",group.df=data.frame(sex=c("M","F")),
covariates="weight")

Comments spanning lines in the .inp file are ignored and deleted as are blank lines. If each line has a unique identifier in the comments then by setting use.comments=TRUE, the text of the comment (e.g.,tag number) will be assigned as the row name in the RMark dataframe. This will only work if each line only represents a single animal or a set of animals in a single group. If file was structured as follows:

 MARK encounter history file (in make believe test.inp): 1001
1 0 10 /*1*/; 1101 0 2 5 /*2*/; 0101 3 1 6 /*3*/; 

an error would occur

 Error in convert.inp("test", group.df = data.frame(sex =
c("M", "F")),: Row names not unique. Set use.comments to default value FALSE

because it would try to use "3" as the row name for the 3 males and the 1 female represented by the last row.

The extension .inp is optional for files with that extension. If the file has a different extension the entire filename must be specified.

Note that there are limitations to this function. You cannot have extra blank lines in the file, the number of columns (tab, space or comma delimited) must be the same in each line unless the line is just a comment line /* */. In the latter case, the /* must begin the line and the */ must end the line with no extra characters (blanks included) in before or after.

Value

Dataframe with fields ch(character encounter history), freq (frequency of encounter history), followed by grouping variables (if any) and then covariates (if any)

Author(s)

Jeff Laake

See Also

process.data

Examples

# This example is excluded from testing to reduce package check time
# MARK example input file
pathtodata=paste(path.package("RMark"),"extdata",sep="/")
dipper=convert.inp(paste(pathtodata,"dipper",sep="/"),
                    group.df=data.frame(sex=c("M","F")))
# Example input files that accompany the MARK electronic book 
#  \url{http://www.phidot.org/software/mark/docs/book/}
bd=convert.inp(paste(pathtodata,"blckduck",sep="/"),
         covariates=c("age","weight","winglen","ci"),use.comments=TRUE)
aa=convert.inp(paste(pathtodata,"aa",sep="/"),
      group.df=data.frame(sex=c("Poor","Good")))
adult=convert.inp(paste(pathtodata,"adult",sep="/"))
age=convert.inp(paste(pathtodata,"age",sep="/"))
age_ya=convert.inp(paste(pathtodata,"age_ya",sep="/"),
      group.df=data.frame(age=c("Young","Adult")))
capsid=convert.inp(paste(pathtodata,"capsid",sep="/"))
clogit_demo=convert.inp(paste(pathtodata,"clogit_demo",sep="/"))
deer=convert.inp(paste(pathtodata,"deer",sep="/"))
ed_males=convert.inp(paste(pathtodata,"ed_males",sep="/"))
F_age=convert.inp(paste(pathtodata,"f_age",sep="/"))
indcov1=convert.inp(paste(pathtodata,"indcov1",sep="/"),
         covariates=c("cov1","cov2"))
indcov2=convert.inp(paste(pathtodata,"indcov2",sep="/"),
          covariates=c("cov1","cov2"))
island=convert.inp(paste(pathtodata,"island",sep="/"))
linear=convert.inp(paste(pathtodata,"linear",sep="/"))
young=convert.inp(paste(pathtodata,"young",sep="/"))
transient=convert.inp(paste(pathtodata,"transient",sep="/"))
ms_gof=convert.inp(paste(pathtodata,"ms_gof",sep="/"))
m_age=convert.inp(paste(pathtodata,"m_age",sep="/"))
ms_cjs=convert.inp(paste(pathtodata,"ms_cjs",sep="/"))
ms_directional=convert.inp(paste(pathtodata,"ms_directional",sep="/"))
ed=convert.inp(paste(pathtodata,"ed",sep="/"),
            group.df=data.frame(sex=c("Male","Female")))
multigroup=convert.inp(paste(pathtodata,"multi_group",sep="/"),
            group.df=data.frame(sex=c(rep("Female",2),rep("Male",2)),
            Colony=rep(c("Good","Poor"),2)))
LD1=convert.inp(paste(pathtodata,"ld1",sep="/"),
           group.df=data.frame(age=c("Young","Adult")))
yngadt=convert.inp(paste(pathtodata,"yngadt",sep="/"),
            group.df=data.frame(age=c("Young","Adult")))
effect_size=convert.inp(paste(pathtodata,"effect_size",sep="/"),
             group.df=data.frame(colony=c("Poor","Good")))
effect_size3=convert.inp(paste(pathtodata,"effect_size3",sep="/"),
             group.df=data.frame(colony=c("1","2","3")))

Description

Computes real parameters from link values

Usage

convert.link.to.real(x, model = NULL, links = NULL, fixed = NULL)

Arguments

Details

Computation of the real parameter from the link value is relatively straightforward for most links and the function inverse.link is used. The only exception is parameters that use the mlogit link which requires the transformation across sets of parameters. This is a convenience function that does the necessary work to convert from link to real for any set of parameters. The appropriate links are obtained from model$links unless the argument links is specified and they will over-ride those in model.

Value

vector of real parameter values

Author(s)

Jeff Laake

See Also

inverse.link,compute.real

Description

Computes real estimates for a dataframe of covariate values and the var-cov matrix of the real estimates.

Usage

covariate.predictions(
  model,
  data = NULL,
  indices = NULL,
  drop = TRUE,
  revised = TRUE,
  mata = FALSE,
  normal.lm = FALSE,
  residual.dfs = 0,
  alpha = 0.025,
  ...
)

Arguments

Details

This function has a similar use as compute.real except that it is specifically designed to compute real parameter estimates for multiple covariate values for either a single model or to compute model averaged estimates across a range of models within a marklist. This is particularly useful for computing and plotting the real parameter as a function of the covariate with pointwise confidence bands (see example below). The function also computes a variance-covariance matrix for the real parameters. For example, assume you had a model with two age classes of young and adult and survial for young was a function of weight and you wanted to estimate survivorship to some adult age as a function of weight. To do that you need the survival for young as a function of weight, the adult survival, the variance of each and their covariance. This function will allow you to accomplish tasks like these and many others.

When a variance-covariance matrix is computed for the real parameters, it can get too large for available memory for a large set of real parameters. Most models contain many possible real parameters to accomodate the general structure even if there are very few unique ones. It is necessary to use the most general structure to accomodate model averaging. Most of the time you will only want to compute the values of a limited set of real parameters but possibly for a range of covariate values. Use the argument indices to select the real parameters to be computed. The index is the value that the real parameter has been assigned with the all-different PIM structure. If you looked at the row numbers in the design data for the dipper example, you would see that the parameter for p and Phi are both numbered 1 to 21. But to deal with multiple parameters effectively they are given a unique number in a specific order. For the CJS model, p follows Phi, so for the dipper example, Phi are numbered 1 to 21 and then p are numbered 22 to 42. You can use the function PIMS to lookup the parameter numbers for a parameter type if you use simplified=FALSE. For example, with the dipper data, you could enter PIMS(dipper.model,"p",simplified=FALSE) and you would see that they are numbered 22 to 42. Alternatively, you can use summary.mark with the argument se=TRUE to see the list of indices named all.diff.index. They are included in a dataframe for each model parameter which enables them to be selected based on the attached data values (e.g., time, group etc). For example, if you fitted a model called dipper.model then you could use summary(dipper.model,se=TRUE)$real to list the indices for each parameter.

The argument data is a dataframe containing values for the covariates used in the models. The names for the fields should match the names of the covariates used in the model. If a time-varying covariate is used then you need to specify the covariate name with the time index included as it is specified in the data. You do not need to specify all covariates that were used. If a covariate in one or more models is not included in data then the mean value will be used for each missing covariate. That can be very useful when you are only interested in prediction for one type of parameters (eg Phi) when there are many covariates that are not interesting in another parameter (e.g., p). For each row in data, each parameter specified in indices is computed with those covariate values. So if there were 5 rows in data and 10 parameters were specified there would be 10 sets of 5 (50) estimates produced. If you do not want the full pairing of data and estimates, create a field called index in data and the estimate for that parameter will be computed with those specific values. For example, if you wanted parameter 1 to be computed with 5 different values of a covariate and then parameter 7 with 2 different covariate values, you could create a dataframe with 5 rows each having an index with the value 1 along with the relevant covariate data and an additional 2 rows with the index with the value 7. If you include the field index in data then you do not need to give a value for the argument indices. However, if you are making the computations for parameters that use an mlogit link you must use the separate indices argument. If you try to use the data.frame(index=...,cov) approach with mlogit parameters and you have covariate values, the function will stop with an error. Also, if you only include a portion of the indices in an mlogit set, it will also stop and issue an error and tell you the set of indices that should be included for that mlogit set. If you were allowed to exclude some indices the result would be incorrect.

Value

A list is returned containing a dataframe of estimates, a var-cov matrix, and a reals list:

Author(s)

Jeff Laake

See Also

compute.real,model.average

Examples

pathtodata=paste(path.package("RMark"),"extdata",sep="/")

# This example is excluded from testing to reduce package check time
#
# indcov1.R
#
# CJS analysis of the individual covariate data from 12.2 of
# Cooch and White
#
# Import data (indcov1.inp) and convert it from the MARK inp file 
# format to the RMark format using the function convert.inp  
# It is defined with 1 group but 2 individual covariates of mass and 
# sqmass
#
indcov1=convert.inp(paste(pathtodata,"indcov1",sep="/"),
               covariates=c("mass","sqmass"))
#
# Next create the processed dataframe and the design data.
#
  ind1.process=process.data(indcov1,model="CJS")
  ind1.ddl=make.design.data(ind1.process)
#
# Next create the function that defines and runs the set of models
# and returns a marklist with the results and a model.table.  
# It does not have any arguments but does use the ind1.process 
# and ind1.ddl objects created above in the workspace. The function 
# create.model.list is the function that creates a dataframe of the
# names of the parameter specifications for each parameter in that 
# type of model. If none are given for any parameter, the default
# specification will be used for that parameter in mark.  The 
# first argument of mark.wrapper is the model list of parameter 
# specifications.  Remaining arguments that are passed to
# mark must be specified using the argument=value specification 
# because the arguments of mark were not repeated in mark.wrapper 
# so they must be passed using the argument=value syntax.
#
ind1.models=function()
{
  Phi.dot=list(formula=~1)
  Phi.mass=list(formula=~mass)
  Phi.mass.plus.mass.squared=list(formula=~mass + sqmass)
  p.dot=list(formula=~1)
  cml=create.model.list("CJS")
  results=mark.wrapper(cml,data=ind1.process,ddl=ind1.ddl,adjust=FALSE,delete=TRUE)
  return(results)
}
#
# Next run the function to create the models and store the results in
# ind1.results which is a marklist. Note that beta estimates will differ
# from Cooch and White results because we are using covariate values 
# directly rather than standardized values.
#
ind1.results=ind1.models()
#
# Next compute real parameter values for survival as a function of 
# mass which are model-averaged over the fitted models.
#
minmass=min(indcov1$mass)
maxmass=max(indcov1$mass)
mass.values=minmass+(0:30)*(maxmass-minmass)/30
Phibymass=covariate.predictions(ind1.results,
   data=data.frame(mass=mass.values,sqmass=mass.values^2),
    indices=c(1))
#
# Plot predicted model averaged estimates by weight with pointwise 
# confidence intervals
#
plot(Phibymass$estimates$mass, Phibymass$estimates$estimate,
   type="l",lwd=2,xlab="Mass(kg)",ylab="Survival",ylim=c(0,.65))
lines(Phibymass$estimates$mass, Phibymass$estimates$lcl,lty=2)
lines(Phibymass$estimates$mass, Phibymass$estimates$ucl,lty=2)

# indcov2.R
#
# CJS analysis of the individual covariate data from 12.3 of 
# Cooch and White
#
# Import data (indcov2.inp) and convert it from the MARK inp file 
# format to the RMark format using the function convert.inp  It is
# defined with 1 group but 2 individual covariates of mass and 
# sqmass
#
indcov2=convert.inp(paste(pathtodata,"indcov2",sep="/"),
        covariates=c("mass","sqmass"))
#
# Standardize covariates
#
actual.mass=indcov2$mass
standardize=function(x,z=NULL)
{
   if(is.null(z))
   {
       return((x-mean(x))/sqrt(var(x)))
   }else 
   {
      return((x-mean(z))/sqrt(var(z)))
   }
}
indcov2$mass=standardize(indcov2$mass)
indcov2$sqmass=standardize(indcov2$sqmass)
#
# Next create the processed dataframe and the design data.
#
  ind2.process=process.data(indcov2,model="CJS")
  ind2.ddl=make.design.data(ind2.process)
#
# Next create the function that defines and runs the set of models and 
# returns a marklist with the results and a model.table.  It does not
# have any arguments but does use the ind1.process and ind1.ddl 
# objects created above in the workspace. The function create.model.list 
# is the function that creates a dataframe of the names of the parameter
# specifications for each parameter in that type of model. If none are 
# given for any parameter, the default specification will be used for
# that parameter in mark.  The first argument of mark.wrapper is the
# model list of parameter specifications.  Remaining arguments that are
# passed to mark must be specified using the argument=value specification 
# because the arguments of mark were not repeated in mark.wrapper so 
# they must be passed using the argument=value syntax.
#
ind2.models=function()
{
  Phi.dot=list(formula=~1)
  Phi.time=list(formula=~time)
  Phi.mass=list(formula=~mass)
  Phi.mass.plus.mass.squared=list(formula=~mass + sqmass)
  Phi.time.x.mass.plus.mass.squared=
        list(formula=~time:mass + time:sqmass)
  Phi.time.mass.plus.mass.squared=
     list(formula=~time*mass + sqmass+ time:sqmass)
  p.dot=list(formula=~1)
  cml=create.model.list("CJS")
  results=mark.wrapper(cml,data=ind2.process,ddl=ind2.ddl,adjust=FALSE,threads=2,delete=TRUE)
  return(results)
}
#
# Next run the function to create the models and store the results in
# ind2.results which is a marklist. Note that beta estimates will differ
# because we are using covariate values directly rather than 
# standardized values.
#
ind2.results=ind2.models()
#
# Next compute real parameter values for survival as a function of 
# mass which are model-averaged over the fitted models. They are 
# standardized individually so the values have to be chosen differently.
#
minmass=min(actual.mass)
maxmass=max(actual.mass)
mass.values=minmass+(0:30)*(maxmass-minmass)/30
sqmass.values=mass.values^2
mass.values=standardize(mass.values,actual.mass)
sqmass.values=standardize(sqmass.values,actual.mass^2)
Phibymass=covariate.predictions(ind2.results,
 data=data.frame(mass=mass.values,sqmass=sqmass.values),
  indices=c(1:7))
#
# Plot predicted model averaged estimates by weight with pointwise 
# confidence intervals
#
par(mfrow=c(4,2))
for (i in 1:7)
{
   mass=minmass+(0:30)*(maxmass-minmass)/30
   x=Phibymass$estimates
   plot(mass,x$estimate[x$par.index==i],type="l",lwd=2,
    xlab="Mass(kg)",ylab="Survival",ylim=c(0,1),main=paste("Time",i))
   lines(mass, x$lcl[x$par.index==i],lty=2)
   lines(mass, x$ucl[x$par.index==i],lty=2)
}

Description

Data and Script to simulate the MSCRD example of 15.7.1 from the MARK book Cooch and White

Format

A data frame with 557 observations on the following 2 variables.

ch

a character vector of encounter histories

freq

a numeric vector of frequencies of each history

Source

This example was constructed by Andrew Paul who is with Fish and Wildlife Division of the Alberta Provincial Government

References

For Cooch and White book see http://www.phidot.org/software/mark/

Examples

# This example is excluded from testing to reduce package check time
#Script to simulate the MSCRD 
#example of 15.7.1 from the MARK
#book
#created by AJP 21 Dec 2010


#convert .inp data - only needed to create crdms
#ch.data<-convert.inp("rd_simple1.inp")
data(crdms)
#set time intervals
#4 primary periods each with 3 secondary occasions
t.int<-c(rep(c(0,0,1),3),c(0,0))

#process data for RMark
crdms.data<-process.data(crdms,model="CRDMS",time.interval=t.int,
			strata.labels=c("1","U"))
#change Psi parameters that are obtained by subtraction
crdms.ddl<-make.design.data(crdms.data,
		parameters=list(Psi=list(subtract.stratum=c("1","1"))))

#create grouping index for unobserved p and c (i.e., always zero)
#up=as.numeric(row.names(crdms.ddl$p[crdms.ddl$p$stratum=="U",]))
# changed 16 Oct 2020 - jll
crdms.ddl$p$fix=ifelse(crdms.ddl$p$stratum=="U",0,NA)
crdms.ddl$c$fix=ifelse(crdms.ddl$c$stratum=="U",0,NA)

#create grouping index to fix Psi for unobs to unbos at time 1
#this isn't necessary but it allows this Psi to be fixed to a value
#that can be flagged and not erroneously interpreted
Psiuu1=as.numeric(row.names(crdms.ddl$Psi[crdms.ddl$Psi$stratum=="U"&
			crdms.ddl$Psi$time==1,]))

#create dummy variable for constraining last Psi in Markovian model
#variable is called ctime for constrained time
crdms.ddl$Psi$ctime=crdms.ddl$Psi$time
crdms.ddl$Psi$ctime[crdms.ddl$Psi$time==3]=2
do_example=function()
{
#Initial assumptions
S.dot=list(formula=~1)  #S equal for both states and constant over time
# p.session=list(formula=~session, share=TRUE,  #p=c varies with session 
# 		fixed=list(index=up,value=0)) #p set to zero for unobs
# changed 16 Oct 2020 - jll
p.session=list(formula=~session, share=TRUE)

#Model 1 - Markovian movement
Psi.markov=list(formula=~ctime+stratum,
		fixed=list(index=Psiuu1,value=9e-99)) #9e-99 is a flag
model.1=mark(crdms.data,ddl=crdms.ddl,
	model.parameters=list(S=S.dot,
		p=p.session,
		Psi=Psi.markov),threads=2,delete=TRUE)

#Model 2 - Random movement
Psi.rand=list(formula=~time)
model.2=mark(crdms.data,ddl=crdms.ddl,
	model.parameters=list(S=S.dot,
		p=p.session,
		Psi=Psi.rand),threads=2,delete=TRUE)

#Model 3 - No movement
Psi.fix=list(formula=~1,fixed=0)
model.3=mark(crdms.data,ddl=crdms.ddl,
	model.parameters=list(S=S.dot,
		p=p.session,
		Psi=Psi.fix),threads=2,delete=TRUE)
		
#collect and store models
crdms.res<-collect.models()

print(crdms.res)
invisible()
}
do_example()

Description

Reads in mcmc file from program MARK in binary and returns an mcmc object that can be used with coda functions which are most easily accessed via codamenu().

Usage

create.mark.mcmc(
  filename,
  ncovs,
  nmeans,
  ndesigns,
  nsigmas,
  nrhos,
  nlogit,
  include = F
)

Arguments

Value

An mcmc object if one chain and an mcmc list if more than one chain. Each can be used with the coda package

Author(s)

Jeff Laake

Description

Creates a dataframe of all combinations of parameter specifications for each parameter in a particular type of MARK model. It is used together with mark.wrapper to run a series of models from sets of parameter specifications.

Usage

create.model.list(model)

Arguments

Details

This function scans the frame of the calling enviroment and collects all list objects that contain a formula and have names that match parameter. where parameter is the name of a type of parameter in the model type. For example, it looks for Phi. and p. for model="CJS". Any number of characters can follow the period. Each of the named objects should specify a list that matches the structure of a parameter specification as described in make.mark.model. It only collects list objects that contain an element named formula, thus it will not collect one like Phi.fixed=list(fixed=1). If you want to do something like that, specify it as Phi.fixed=list(formula=~1,fixed=1). It is safest to use this inside a function that defines all of the parameter specifications as shown in the example below. The primary use for this function is to create a dataframe which is passed to mark.wrapper to construct and run each of the models. It was written as a separate function to provide flexibility to add/delete/modify the list prior to passing to mark.wrapper. For example, only certain combinations may make sense for some parameter specifications. Thus you could define a set to create all the combinations and then delete the ones from the dataframe that do not make sense. you want, add others and re-run the function and merge the resulting dataframes. If there are no specifications found for a particular model parameter, it is not included in the list and when it is passed to make.mark.model, the default specification will be used for that parameter.

Value

dataframe of all combinations of parameter specifications for a model. Each field (column) is the name of a type of parameter (e.g., p and Phi for CJS). The values are character strings identifying particular parameter specifications.

Author(s)

Jeff Laake

See Also

mark.wrapper

Examples

# This example is excluded from testing to reduce package check time
#
# Compare this to the run.dipper shown under ?dipper
# It is only necessary to create each parameter specification and 
# create.model.list and mark.wrapper will create and run models for 
# each combination. Notice that the naming of the parameter 
# specifications has been changed to accommodate format for 
# create.model.list. Only a subset of the parameter specifications 
# are used here in comparison to other run.dipper
#
data(dipper)
run.dipper=function()
{
     #
     # Process data
     #
     dipper.processed=process.data(dipper,groups=("sex"))
     #
     # Create default design data
     #
     dipper.ddl=make.design.data(dipper.processed)
     #
     # Add Flood covariates for Phi and p that have different values
     #
     dipper.ddl$Phi$Flood=0
     dipper.ddl$Phi$Flood[dipper.ddl$Phi$time==2 |
                  dipper.ddl$Phi$time==3]=1
     dipper.ddl$p$Flood=0
     dipper.ddl$p$Flood[dipper.ddl$p$time==3]=1
     #
     #  Define range of models for Phi
     #
     Phi.dot=list(formula=~1)
     Phi.time=list(formula=~time)
     Phi.sex=list(formula=~sex)
     Phi.Flood=list(formula=~Flood)
     #
     #  Define range of models for p
     #
     p.dot=list(formula=~1)
     p.time=list(formula=~time)
     p.Flood=list(formula=~Flood)
     #
     # Run all pairings of models
     #
     dipper.model.list=create.model.list("CJS")
     dipper.results=mark.wrapper(dipper.model.list,
              data=dipper.processed,ddl=dipper.ddl,delete=TRUE)
     #
     # Return model table and list of models
     #
     return(dipper.results)
}
dipper.results=run.dipper()

Description

This data represents a set of independent double observer road-transect survey data of white-tailed deer on Brosnan Forest, South Carolina surveyed in August, 2005-2009. The primary reason for this package is to provide a completely reproducible example of the analysis from Collier et al. (2012). We used a Huggins closed capture model implemented in MARK http://www.phidot.org/software/mark/ via RMark both of which will need to be installed on the system to use this package. The data have 2 time periods (primary observer (t1) was a thermal imager, secondary observer (t2) was a spotlight observer in the same vehicle on the same side) with the primary objective of the study being to evaluate the detection (recapture) rates of white-tailed deer using spotlights as a survey method.

Format

The format is a data frame with 4508 observations on the following 7 variables.

SL (spotlight)

0/1 whether deer was missed/seen by the spotlight observer

TI (thermal imager)

0/1 whether deer was missed/seen by the thermal imager observer

Group

Factor with 79 levels representing each unique paired (TI-SL) survey conducted

Year

Factor with 5 levels for year of survey

MaxCount

Count of maximum number of deer seen for each survey, only needed for bootstrapp analysis in MARK, not used in bfdeeR package

Cluster

Value assigning each deer to a specific observation cluster, only needed for bootstrapp analysis in MARK, not used in bfdeeR package

MgmtUnit

Management unit identification

Details

In addition to detailing the analysis used by Collier et al. (2012), this example documents the use of the share argument in the RMark parameter specification because there is presently very little documentation on the use of share. Parameters in MARK models rarely share columns of the design matrix. For example while you might want to use the same covariate for survival and capture probability, you would never use the same beta (same column of the design matrix) for each parameter. However, there are exceptions when the parameters represent similar quantities and that is when the share argument is useful. For example, in the closed capture models p is initial capture probability and c is recapture probability. In this case, it would make perfect sense to use the same column of the design matrix for both parameters. The most obvious case is to fit a model in which p=c.

In RMark, certain pairs of parameters have been identified as similar and shareable. These can be found in the file parameters.txt which is in the RMark directory in your R library. With each pair that is shareable, the first one listed is the primary parameter. When you want to share columns in the design matrix, share=TRUE is added to the specification of the primary parameter. A parameter specification is not given for the other secondary parameter when they are shared. When RMark, sees that the parameters are to be shared it creates a pooled set of design data and adds a column with the name of the secondary parameter and its value is 0 for the rows for the primary parameter and 1 for the rows for the secondary parameter. For example, with the closed capture model if share=TRUE is added to the parameter specification for p, a model is not specified for c, and the pooled design data set contains a field called c. The added field allows construction of models where there are restricted differences between the parameters. For example, p=list(formula=~time+c,share=TRUE) will fit a model in which capture probability varies by time and recapture probability includes an additive difference on the link scale. Because the design data are pooled when you share parameters, if you modify design data for one of the parameters, the other most be modified as as well, so the columns of the design data for both parameters are the same or RMark will give an error.

The argument share is used in all the candidate models in the below example analysis. As a simplified example of how share works, look at the candidate models in the bfrun{} function call named mod.2 and mod.2a (note that mod.2a was not included in the supplemental file available from the Journal of Wildlife Management and is only included in this package). Both of these models are conducting the exact same analysis, with the first mod.2, we used the formula ~time (if you don't know what this means go read the MARKBOOK at http://www.phidot.org/software/mark/. Notice, however, we used the argument share in mod.2, which tells RMark to share columns of the MARK design matrix. For comparison, so you can evaluate how share works for yourself, mod.2a recreates the same analysis as mod.2, but uses the approach more typical to MARK analyses where each parameter is specified independently and uniquely.

Author(s)

Bret Collier

References

Collier, B. A., S. S. Ditchkoff, J. B. Raglin, and C. R. Ruth. 2012. Spotlight surveys for white-tailed deer: monitoring panacea or exercise in futility? Journal of Wildlife Management, In Press.

Examples

# This example is excluded from testing to reduce package check time
data(deer)
x=data.frame(ch=paste(deer$TI, deer$SL, sep=""), Survey=factor(deer$Group), 
     Year=factor(deer$Year), Cluster=deer$Cluster, MgtUnit=factor(deer$MgmtUnit))
x$ch=as.character(x$ch)
bfrun=function(){
x.proc=process.data(x, model="Huggins", groups=c("Survey", "Year", "MgtUnit"))
x.ddl=make.design.data(x.proc)

#Silly Null model, constant p & c sharing 1 parameter (one detection estimate)
p.shared=list(formula=~1,share=TRUE)
mod.1=mark(x.proc, x.ddl, model.parameters=list(p=p.shared), invisible=FALSE,delete=TRUE)
 
#2 Parameter Null Model, constant p, constant c, different p and c (one estimate for each; p ne c)
#p(time), c(-), share=TRUE, detection is time dependent, with recapture parameter shared
p.sharetime=list(formula=~time, share=TRUE)
mod.2=mark(x.proc, x.ddl, model.parameters=list(p=p.sharetime), invisible=FALSE,delete=TRUE)

#2a Parameter Null Model, constant p, constant c,
# different p and c (one estimate for each; p ne c) not using share
mod.2a=mark(x.proc, x.ddl, model.parameters=list(p=list(formula=~1), c=list(formula=~1)),
            delete=TRUE)

#Fully parameterized model, different p and c for each survey transect replicate, 
# management unit, method (TI or SL) and any observers
p.survey=list(formula=~Survey*time, share=TRUE)
mod.3=mark(x.proc, x.ddl, model.parameters=list(p=p.survey), invisible=FALSE,delete=TRUE)

#p(MU), c(MU), initial detection and recapture differ and are management unit dependent
p.mu=list(formula=~MgtUnit*time, share=TRUE)
mod.4=mark(x.proc, x.ddl, model.parameters=list(p=p.mu), invisible=FALSE,delete=TRUE)

#p(MU) detection is management unit dependent
p.mu=list(formula=~MgtUnit, share=TRUE)
mod.5=mark(x.proc, x.ddl, model.parameters=list(p=p.mu), invisible=FALSE,delete=TRUE)

#p(Yr + MgtUnit),  detection is year + MgtUnit
p.yearMgtUnit=list(formula=~Year*time+MgtUnit, share=TRUE)
mod.6=mark(x.proc, x.ddl, model.parameters=list(p=p.yearMgtUnit), invisible=FALSE,delete=TRUE)

#p(Year), initial detection and recapture are year dependent
p.year=list(formula=~Year*time, share=TRUE)
mod.7=mark(x.proc, x.ddl, model.parameters=list(p=p.year), invisible=FALSE,delete=TRUE)

return(collect.models())
}
bf.out=bfrun()
bf.out

#export function to send dataset and covariates data to MARK for bootstrap analysis 
#(not run but here for completeness)
#export.MARK(x.proc, "BFdeer", mod.3, replace=TRUE, ind.covariates="all")

Description

This function computes the delta method std errors or v-c matrix for a sum, cumsum (vector of cummulative sums), prod (product), or cumprod(vector of cummulative products) of a set of estimates.

Usage

deltamethod.special(function.name, mean, cov, ses = TRUE)

Arguments

Details

This function computes the delta method std errors or v-c matrix for a sum, cumsum (vector of cummulative sums), prod (product), or cumprod(vector of cummulative products). It uses the function deltamethod from the msm package and constructs the necessary formula for these special cases. It will load the msm pacakge but assumes that it has already been installed. See the msm documentation for a complete description on how the deltamethod function works. If ses=TRUE, it returns a vector of std errors for each of functions of the estimates contained in mean. If ses=F, then it returns a v-c matrix for the functions of the estimates contained in mean. cov is the input v-c matrix of the estimates.

Value

either a vector of standard errors (ses=TRUE) or a variance-covariance matrix (ses=FALSE)

Author(s)

Jeff Laake

Examples

# This example is excluded from testing to reduce package check time
#
# The following are examples only to demonstrate selecting different 
# model sets for adjusting chat and showing model selection table. 
# It is not a realistic analysis.
#
  data(dipper)
  mod=mark(dipper,model.parameters=list(Phi=list(formula=~time)),delete=TRUE)
  rr=get.real(mod,"Phi",se=TRUE,vcv=TRUE)
  deltamethod.special("prod",rr$estimates$estimate[1:6],rr$vcv.real)
  deltamethod.special("cumprod",rr$estimates$estimate[1:6],rr$vcv.real,ses=FALSE)
  deltamethod.special("sum",rr$estimates$estimate[1:6],rr$vcv.real)
  deltamethod.special("cumsum",rr$estimates$estimate[1:6],rr$vcv.real,ses=FALSE)

Description

The annotated code below is a companion to A Gentle Introduction to Program MARK, Chapter 20: Density Estimation (http://www.phidot.org/software/mark/docs/book/pdf/chap20.pdf). It requires the file "Density.txt", which is the RMark analog to the example Density Estimation input file distributed with MARK. These are simulated data intended to mimic a study of small mammals, such as deer mice, sampled at 2 sites (habitat types), A and B. Each habitat type was sampled with a (10 x 10) live-trapping grid (10m trap spacing). There are 5 occasions. In addition to marking each mouse with an individually identifiable ear tag, 50 percent of the individuals captured were fitted with a small VHF transmitter. These radio-tagged individuals were located once during the day and once at night for 5 days immediately after mark-recapture sampling (n = 10 locations total per animal) and each location was recorded as in or out of the study site. The single covariate we recorded is the distance to the edge (DTE) of the of the site from the mean trap location of each individual (i.e., compute the mean trap location for each individual captured >1 time, then compute the minimum distance from this mean location to the edge of the site).

Format

A data frame with 32 observations of 5 variables

ch

a character vector containing the capture history for 5 occasions

TotalLocations

The total number of telemetry locations if telemetered; otherwise a .

TotalIn

total number of locations in the original site if telemetered; otherwise a .

Site

the original site A or B

DTE

distance to the edge of the site when originally caught

Author(s)

Jake Ivan

Examples

#Read in Density Estimation input file specific to RMark

#Add 2 covariates that will be used for threshhold models - see below & p. 20-14, 20-15
#Specify data type - use "Densitypc" for this example, which is "Density Estimation with
# Huggins p and c". Could also use "DensityRanpc" (Huggins p and c with random effects),
#"DensityHet" (Huggins heterogeniety with pi and p), "DensityFHet" (Huggins full
#heterogeneity with pi and p) and DensityFHet (Huggins Full heterogeniety with pi, p, and c).
#Be sure to specify areas argument in process.data for this model. It will not run if you don't
# give it the area of each study site

data(Density)
#Create variables for threshhold model - see below & p. 20-14, 20-15
Density$Thresh15 <- ifelse(Density$DTE<15, Density$DTE, 15)
#Create variables for threshhold model - see below & p. 20-14, 20-15
Density$Thresh25 <- ifelse(Density$DTE<25, Density$DTE, 25)
data_proc <- process.data(Density, model="Densitypc", groups = c("Site"), areas=rep(0.81,2))
data_ddl <- make.design.data(data_proc)

#Run model p(.)p~(.) from p. 20-9, 20-10. View results.
p_dot <- list(formula = ~1, share=TRUE)
ptilde_dot <- list(formula = ~1)
model1 <- mark(data_proc,data_ddl,model.parameters=list(p=p_dot,ptilde=ptilde_dot),delete=TRUE)

#Run models p(site)p~(.) and p(.)p~(site) as indicated on p. 20-11. View results.
p_site <- list(formula = ~1 + Site, share=TRUE)
ptilde_site <- list(formula = ~1 + Site)
model2 <- mark(data_proc, data_ddl, model.parameters=list(p=p_dot, ptilde=ptilde_site),
                  delete=TRUE)
model3 <- mark(data_proc, data_ddl, model.parameters=list(p=p_site, ptilde=ptilde_dot),
                    delete=TRUE)


#Run model p(DTE)p~(DTE) as indicated on p. 20-12. View results.
p_DTE <- list(formula = ~1 + DTE, share=TRUE)
ptilde_DTE <- list(formula= ~1 + DTE)
model4 <- mark(data_proc, data_ddl, model.parameters = list(p=p_DTE, ptilde=ptilde_DTE),
                  delete=TRUE)


#Compute Model Selection Table that appears on p. 20-12. View results.
ModSelTable <- collect.models(type="Densitypc")
ModSelTable

#Run Threshhold models from p. 20-15.
p_DTE_Thresh15 <- list(formula = ~1 + Thresh15, share=TRUE)
p_DTE_Thresh25 <- list(formula = ~1 + Thresh25, share=TRUE)
model5 <- mark(data_proc, data_ddl, model.parameters = list(p=p_DTE_Thresh15, ptilde=ptilde_DTE)
               ,delete=TRUE)
model6 <- mark(data_proc, data_ddl, model.parameters = list(p=p_DTE_Thresh25, ptilde=ptilde_DTE)
               ,delete=TRUE)


#Re-compute Model Selection Table that appears on p. 20-16
ModSelTable <- collect.models(type="Densitypc")
ModSelTable

Description

Computes derivatives of inverse of link functions (real estimates) with respect to the beta parameters which are used for delta method computation of the var-cov matrix of real parameters.

Usage

deriv_inverse.link(real, x, link)

Arguments

Details

Note: that function was renamed to deriv_inverse.link to avoid S3 generic class conflicts. The derivatives of the inverse of the link functions are simple computations using the real values and the design matrix values. The body of the function is as follows:

switch(link, logit=x*real*(1-real), log=x*real,
loglog=-real*x*log(real), cloglog=-log(1-real)*x*(1-real), identity=x,
mlogit=x*real*(1-real)) 

Value

Vector of derivative values computed at values of real parameters

Author(s)

Jeff Laake

See Also

inverse.link, compute.real

Description

A capture-recapture data set on European dippers from France that accompanies MARK as an example analysis using the CJS and POPAN models. The dipper data set was orginally described as an example by Lebreton et al (1992).

Format

A data frame with 294 observations on the following 2 variables.

ch

a character vector containing the encounter history of each bird

sex

the sex of the bird: a factor with levels Female Male

Details

This is a data set that accompanies program MARK as an example for CJS and POPAN analyses. The data can be stratified using sex as a grouping variable. The functions run.dipper, run.dipper.alternate, run.dipper.popan defined below in the examples mimic the models used in the dbf file that accompanies MARK. Note that the models used in the MARK example use PIM coding with the sin link function which is often better at identifying the number of estimable parameters. The approach used in the R code uses design matrices and cannot use the sin link and is less capable at counting parameters. These differences are illustrated by comparing the results of run.dipper and run.dipper.alternate which fit the same set of "CJS" models. The latter fits the models with constraints on some parameters to achieve identifiability and the former does not. Although it does not influence the selection of the best model it does infleunce parameter counts and AIC ordering of some of the less competitive models. In using design matrices it is best to constrain parameters that are confounded (e.g., last occasion parameters in Phi(t)p(t) CJS model) when possible to achieve more reliable counts of the number of estimable parameters. See adjust.parameter.count for more dicussion on this point.

Note that the covariate "sex" defined in dipper has values "Male" and "Female". It cannot be used directly in a formula for MARK without using it do define groups because MARK.EXE will be unable to read in a covariate with non-numeric values. By using groups="sex" in the call the process.data a factor "sex" field is created that can be used in the formula. Alternatively, a new covariate could be defined in the data with say values 0 for Female and 1 for Male and this could be used without defining groups because it is numeric. This can be done easily by translating the values of the coded variables to a numeric variable. Factor variables are numbered 1..k for k levels in alphabetic order. Since Female < Male in alphabetic order then it is level 1 and Male is level 2. So the following will create a numeric sex covariate.

 dipper$numeric.sex=as.numeric(dipper$sex)-1 

See export.chdata for an example that creates a .inp file for MARK with sex being used to describe groups and a numeric sex covariate.

Source

Lebreton, J.-D., K. P. Burnham, J. Clobert, and D. R. Anderson. 1992. Modeling survival and testing biological hypotheses using marked animals: case studies and recent advances. Ecol. Monogr. 62:67-118.

Examples

# This example is excluded from testing to reduce package check time
data(dipper)
dipper.model=mark(dipper,delete=TRUE)
run.dipper=function()
{
#
# Process data
#
dipper.processed=process.data(dipper,groups=("sex"))
#
# Create default design data
#
dipper.ddl=make.design.data(dipper.processed)
#
# Add Flood covariates for Phi and p that have different values
#
dipper.ddl$Phi$Flood=0
dipper.ddl$Phi$Flood[dipper.ddl$Phi$time==2 | dipper.ddl$Phi$time==3]=1
dipper.ddl$p$Flood=0
dipper.ddl$p$Flood[dipper.ddl$p$time==3]=1
#
#  Define range of models for Phi
#
Phidot=list(formula=~1)
Phitime=list(formula=~time)
Phisex=list(formula=~sex)
Phisextime=list(formula=~sex+time)
Phisex.time=list(formula=~sex*time)
PhiFlood=list(formula=~Flood)
#
#  Define range of models for p
#
pdot=list(formula=~1)
ptime=list(formula=~time)
psex=list(formula=~sex)
psextime=list(formula=~sex+time)
psex.time=list(formula=~sex*time)
pFlood=list(formula=~Flood)
#
# Run assortment of models
#
dipper.phidot.pdot          =mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phidot,p=pdot),delete=TRUE)
dipper.phidot.pFlood      	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phidot,p=pFlood),delete=TRUE)
dipper.phidot.psex        	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phidot,p=psex),delete=TRUE)
dipper.phidot.ptime       	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phidot,p=ptime),delete=TRUE)
dipper.phidot.psex.time		=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phidot,p=psex.time),delete=TRUE)
dipper.phitime.ptime      	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phitime, p=ptime),delete=TRUE)
dipper.phitime.pdot       	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phitime,p=pdot),delete=TRUE)
dipper.phitime.psex		=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phitime,p=psex),delete=TRUE)
dipper.phitime.psex.time	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phitime,p=psex.time),delete=TRUE)
dipper.phiFlood.pFlood    	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=PhiFlood, p=pFlood),delete=TRUE)
dipper.phisex.pdot        	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phisex,p=pdot),delete=TRUE)
dipper.phisex.psex        	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phisex,p=psex),delete=TRUE)
dipper.phisex.psex.time        	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phisex,p=psex.time),delete=TRUE)
dipper.phisex.ptime       	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phisex,p=ptime),delete=TRUE)
dipper.phisextime.psextime	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phisextime,p=psextime),delete=TRUE)
dipper.phisex.time.psex.time	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phisex.time,p=psex.time),delete=TRUE)
dipper.phisex.time.psex 	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phisex.time,p=psex),delete=TRUE)
dipper.phisex.time.pdot		=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phisex.time,p=pdot),delete=TRUE)
dipper.phisex.time.ptime	=mark(dipper.processed,dipper.ddl,
                 model.parameters=list(Phi=Phisex.time,p=ptime),delete=TRUE)
#
# Return model table and list of models
#
return(collect.models() )
}

dipper.results=run.dipper()

run.dipper.alternate=function()
{
#
# Process data
#
dipper.processed=process.data(dipper,groups=("sex"))
#
# Create default design data
#
dipper.ddl=make.design.data(dipper.processed)
#
# Add Flood covariates for Phi and p that have different values
#
dipper.ddl$Phi$Flood=0
dipper.ddl$Phi$Flood[dipper.ddl$Phi$time==2 | dipper.ddl$Phi$time==3]=1
dipper.ddl$p$Flood=0
dipper.ddl$p$Flood[dipper.ddl$p$time==3]=1
#
#  Define range of models for Phi
#
Phidot=list(formula=~1)
Phitime=list(formula=~time)
Phitimec=list(formula=~time,fixed=list(time=6,value=1))
Phisex=list(formula=~sex)
Phisextime=list(formula=~sex+time)
Phisex.time=list(formula=~sex*time)
PhiFlood=list(formula=~Flood)
#
#  Define range of models for p
#
pdot=list(formula=~1)
ptime=list(formula=~time)
ptimec=list(formula=~time,fixed=list(time=7,value=1))
psex=list(formula=~sex)
psextime=list(formula=~sex+time)
psex.time=list(formula=~sex*time)
psex.timec=list(formula=~sex*time,fixed=list(time=7,value=1))
pFlood=list(formula=~Flood)
#
# Run assortment of models
#
dipper.phidot.pdot          =mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phidot,p=pdot),delete=TRUE)
dipper.phidot.pFlood      	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phidot,p=pFlood),delete=TRUE)
dipper.phidot.psex        	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phidot,p=psex),delete=TRUE)
dipper.phidot.ptime       	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phidot,p=ptime),delete=TRUE)
dipper.phidot.psex.time		=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phidot,p=psex.time),delete=TRUE)
dipper.phitime.ptimec      	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phitime, p=ptimec),delete=TRUE)
dipper.phitime.pdot       	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phitime,p=pdot),delete=TRUE)
dipper.phitime.psex		=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phitime,p=psex),delete=TRUE)
dipper.phitimec.psex.time	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phitimec,p=psex.time),delete=TRUE)
dipper.phiFlood.pFlood    	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=PhiFlood, p=pFlood),delete=TRUE)
dipper.phisex.pdot        	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phisex,p=pdot),delete=TRUE)
dipper.phisex.psex        	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phisex,p=psex),delete=TRUE)
dipper.phisex.psex.time        	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phisex,p=psex.time),delete=TRUE)
dipper.phisex.ptime       	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phisex,p=ptime),delete=TRUE)
dipper.phisextime.psextime	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phisextime,p=psextime),adjust=FALSE,delete=TRUE)
dipper.phisex.time.psex.timec	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phisex.time,p=psex.timec),delete=TRUE)
dipper.phisex.time.psex 	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phisex.time,p=psex),delete=TRUE)
dipper.phisex.time.pdot		=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phisex.time,p=pdot),delete=TRUE)
dipper.phisex.time.ptimec	=mark(dipper.processed,dipper.ddl,
                  model.parameters=list(Phi=Phisex.time,p=ptimec),delete=TRUE)
#
# Return model table and list of models
#
return(collect.models() )
}
dipper.results.alternate=run.dipper.alternate()
#
# Merge two sets of models into a single model list and include the 
# initial model as a demo for merge.mark
#
dipper.cjs=merge.mark(dipper.results,dipper.results.alternate,dipper.model)
dipper.cjs
#
# next delete some of the models to show how this is done with remove.mark
#
dipper.cjs=remove.mark(dipper.cjs,c(2,4,9))
dipper.cjs

run.dipper.popan=function()
{
#
# Process data
#
dipper.processed=process.data(dipper,model="POPAN",group="sex")
#
# Create default design data
#
dipper.ddl=make.design.data(dipper.processed)
#
# Add Flood covariates for Phi and p that have different values
#
dipper.ddl$Phi$Flood=0
dipper.ddl$Phi$Flood[dipper.ddl$Phi$time==2 | dipper.ddl$Phi$time==3]=1
dipper.ddl$p$Flood=0
dipper.ddl$p$Flood[dipper.ddl$p$time==3]=1
#
#  Define range of models for Phi
#
Phidot=list(formula=~1)
Phitime=list(formula=~time)
Phisex=list(formula=~sex)
Phisextime=list(formula=~sex+time)
Phisex.time=list(formula=~sex*time)
PhiFlood=list(formula=~Flood)
#
#  Define range of models for p
#
pdot=list(formula=~1)
ptime=list(formula=~time)
psex=list(formula=~sex)
psextime=list(formula=~sex+time)
psex.time=list(formula=~sex*time)
pFlood=list(formula=~Flood)
#
#  Define range of models for pent
#
pentsex.time=list(formula=~sex*time)
#
#  Define range of models for N
#
Nsex=list(formula=~sex)
#
# Run assortment of models
#
dipper.phisex.time.psex.time.pentsex.time=mark(dipper.processed,dipper.ddl,
model.parameters=list(Phi=Phisex.time,p=psex.time,pent=pentsex.time,N=Nsex),
invisible=FALSE,adjust=FALSE,delete=TRUE)
dipper.phisex.time.psex.pentsex.time=mark(dipper.processed,dipper.ddl,
model.parameters=list(Phi=Phisex.time,p=psex,pent=pentsex.time,N=Nsex),
invisible=FALSE,adjust=FALSE,delete=TRUE)
#
# Return model table and list of models
#
return(collect.models() )
}

dipper.popan.results=run.dipper.popan()

# *****************************************************************
# Here is an example of user specified links for each real parameter
 data(dipper)
 dipper.proc=process.data(dipper)
 dipper.ddl=make.design.data(dipper.proc)
# dummy run of make.mark.model to get links and design data. 
# parm.specific set to TRUE so it will create a link for 
# each parameter because for this model they are all the
# same (logit) and if this was not specified you'ld get a vector with one element
 dummy=make.mark.model(dipper.proc,dipper.ddl,simplify=FALSE,parm.specific=TRUE)
 input.links=dummy$links
# get model indices for p where time=4
 log.indices=dipper.ddl$p$model.index[dipper.ddl$p$time==4]
# assign those links to log
 input.links[log.indices]="Log"
# Now these can be used with any call to mark
 mymodel=mark(dipper.proc,dipper.ddl,input.links=input.links,delete=TRUE)
 summary(mymodel)

Description

An example occupancy data set used as exercise 7 in the occupancy website developed by Donovan and Hines.

Format

A data frame with 20 observations (sites) on the following 2 variables.

ch

a character vector containing the presence (1) and absence (0) for each of 5 visits to the site

freq

frequency of sites (always 1)

Details

This is a data set from exercise 7 of Donovan and Hines occupancy web site (https://blog.uvm.edu/tdonovan-vtcfwru/occupancy-estimation-and-modeling-ebook/).

Examples

# Donovan.7 can be created with
# Donovan.7=convert.inp("Donovan.7.inp")

do.exercise.7=function()
{
  data(Donovan.7)
# Estimates from following agree with estimates on website but the
# log-likelihood values do not agree.  Maybe a difference in whether the
# constant binomial coefficients are included.
  Donovan.7.poisson=mark(Donovan.7,model="OccupRNPoisson",invisible=FALSE,threads=1,delete=TRUE)
# THe following model was not in exercise 7.
  Donovan.7.negbin=mark(Donovan.7,model="OccupRNNegBin",invisible=FALSE,threads=1,delete=TRUE)
  return(collect.models())
}
exercise.7=do.exercise.7()
# Remove # to see output
# print(exercise.7)

Description

An example occupancy data set used as exercise 8 in the occupancy website developed by Donovan and Hines.

Format

A data frame with 20 observations (sites) on the following 2 variables.

ch

a character vector containing the counts for each of 5 visits to the site

freq

frequency of sites (always 1)

Details

This is a data set from exercise 8 of Donovan and Hines occupancy web site (https://blog.uvm.edu/tdonovan-vtcfwru/occupancy-estimation-and-modeling-ebook/). In MARK, it uses 2 digits to allow a count of 0 to 99 at each site, so the history has 10 digits for 5 visits (occasions).

Examples

# This example is excluded from testing to reduce package check time
# Donovan.8 can be created with
# Donovan.8=convert.inp("Donovan.8.inp")
do.exercise.8=function()
{
  data(Donovan.8)
# Results agree with the values on the website.
  Donovan.8.poisson=mark(Donovan.8,model="OccupRPoisson",invisible=FALSE,threads=2,delete=TRUE)
# The following model was not in exercise 8. The NegBin model does 
# better if it is initialized with the r and lambda from the poisson.
  Donovan.8.negbin=mark(Donovan.8,model="OccupRNegBin",
    initial=Donovan.8.poisson,invisible=FALSE,threads=2,delete=TRUE)
  return(collect.models())
}
exercise.8=do.exercise.8()
# Remove # to see output
# print(exercise.8)

Description

A capture-recapture data set on rabbits derived from Edwards and Eberhardt (1967) that accompanies MARK as an example analysis using the closed population models.

Format

A data frame with 76 observations on the following variable.

ch

a character vector

Details

This data set is used in MARK to illustrate the various closed population models including "Closed", "HetClosed", "FullHet","Huggins","HugHet", and "FullHugHet". The first 3 include N in the likelihood whereas the last 3 are based on the Huggins approach which does not use N in the likelihood. The Het... and FullHet... models are based on the Pledger mixture model approach. Some of the examples demonstrate the use of the share argument in the model.parameters list for parameter p which allows sharing common values for p and c.

Source

Edwards, W.R. and L.L. Eberhardt 1967. Estimating cottontail abundance from live trapping data. J. Wildl. Manage. 31:87-96.

Examples

# This example is excluded from testing to reduce package check time
#
# get data
#
data(edwards.eberhardt)
#
# create function that defines and runs the analyses as defined in 
# MARK example dbf file
#
run.edwards.eberhardt=function()
{
#
#  Define parameter models
#
pdotshared=list(formula=~1,share=TRUE)
ptimeshared=list(formula=~time,share=TRUE)
ptime.c=list(formula=~time+c,share=TRUE)
ptimemixtureshared=list(formula=~time+mixture,share=TRUE)
pmixture=list(formula=~mixture)
#
# Run assortment of models
#
#
#   Capture Closed models
#
#  constant p=c
ee.closed.m0=mark(edwards.eberhardt,model="Closed",
                   model.parameters=list(p=pdotshared),delete=TRUE)
#  constant p and constant c but different
ee.closed.m0c=mark(edwards.eberhardt,model="Closed",delete=TRUE)
#  time varying p=c
ee.closed.mt=mark(edwards.eberhardt,model="Closed",
                   model.parameters=list(p=ptimeshared),delete=TRUE)
#
#  Closed heterogeneity models
#
#  2 mixtures Mh2
ee.closed.Mh2=mark(edwards.eberhardt,model="HetClosed",
                   model.parameters=list(p=pmixture),delete=TRUE)
#  Closed Mth2 - p different for time; mixture additive
ee.closed.Mth2.additive=mark(edwards.eberhardt,model="FullHet",
                   model.parameters=list(p=ptimemixtureshared),adjust=TRUE,delete=TRUE)
#
#    Huggins models
#
# p=c constant over time
ee.huggins.m0=mark(edwards.eberhardt,model="Huggins",
                   model.parameters=list(p=pdotshared),delete=TRUE)
# p constant c constant but different; this is default model for Huggins
ee.huggins.m0.c=mark(edwards.eberhardt,model="Huggins",delete=TRUE)
# Huggins Mt
ee.huggins.Mt=mark(edwards.eberhardt,model="Huggins",
                   model.parameters=list(p=ptimeshared),adjust=TRUE,delete=TRUE)
#
#    Huggins heterogeneity models
#
#  Mh2 - p different for mixture
ee.huggins.Mh2=mark(edwards.eberhardt,model="HugHet",
                   model.parameters=list(p=pmixture),delete=TRUE)
#  Huggins Mth2 - p different for time; mixture additive
ee.huggins.Mth2.additive=mark(edwards.eberhardt,model="HugFullHet",
                   model.parameters=list(p=ptimemixtureshared),adjust=TRUE,delete=TRUE)
#
# Return model table and list of models
#
return(collect.models() )
}
#
# fit models in mark by calling function created above
#
ee.results=run.edwards.eberhardt()

Description

A simulated data set from CJS model to demonstrate use of grouping variables and individual covariates in an analysis of a mark model. The true model for the simulated data is S(age+year)p(year).

Format

A data frame with 6000 observations on the following 5 variables.

ch

a character vector

weight

a numeric vector

age

a factor with levels 1 2 3

sex

a factor with levels F M

region

a factor with levels 1 2 3 4

Details

The weight, age and region are static variables that are defined based on the values when the animal was released. age is a factor variable representing an age level. The actual ages (at time of release) are 0,1,2 for the 3 levels respectively.

Examples

# This example is excluded from testing to reduce package check time
data(example.data)
run.example=function()
{
PhiTime=list(formula=~time)
pTimec=list(formula=~time,fixed=list(time=7,value=1))
pTime=list(formula=~time)
PhiAge=list(formula=~age)
Phidot=list(formula=~1)
PhiweightTime=list(formula=~weight+time)
PhiTimeAge=list(formula=~time+age)
mod1=mark(example.data,groups=c("sex","age","region"),
                           initial.ages=c(0,1,2),delete=TRUE)
mod2=mark(example.data,model.parameters=list(p=pTimec,Phi=PhiTime),
          groups=c("sex","age","region"),initial.ages=c(0,1,2),delete=TRUE)
mod3=mark(example.data,model.parameters=list(Phi=Phidot,p=pTime),
          groups=c("sex","age","region"),initial.ages=c(0,1,2),delete=TRUE)
mod4=mark(example.data,model.parameters=list(Phi=PhiTime),
          groups=c("sex","age","region"),initial.ages=c(0,1,2),delete=TRUE)
mod5=mark(example.data,model.parameters=list(Phi=PhiTimeAge),
          groups=c("sex","age","region"),initial.ages=c(0,1,2),delete=TRUE)
mod6=mark(example.data,model.parameters=list(Phi=PhiAge,p=pTime),
          groups=c("sex","age","region"),initial.ages=c(0,1,2),delete=TRUE)
mod7=mark(example.data,model.parameters=list(p=pTime,Phi=PhiweightTime),
          groups=c("sex","age","region"),initial.ages=c(0,1,2),delete=TRUE)
mod8=mark(example.data,model.parameters=list(Phi=PhiTimeAge,p=pTime),
          groups=c("sex","age","region"),initial.ages=c(0,1,2),delete=TRUE)
return(collect.models())
}
example.results=run.example()

Description

Creates a MARK .inp file from processed data list that can be used to create a MARK .dbf file for use with MARK directly rather than with the RMark package.

Usage

export.chdata(data, filename, covariates = NULL, replace = FALSE)

Arguments

Details

The default is to include none of the covariates in the processed data list. All of the covariates can be passed by setting covariates="all".

After you have created the MARK .dbf/.fpt files with the exported .inp file, then you can use export.chdata to export models that can be imported into the MARK interface. However note the following: ***Warning*** Make sure that you use the .inp created by export.chdata with your processed data to create the MARK .dbf file rather than using a separate similar .inp file. It is essential that the group structure and ordering of groups matches between the .inp file and the exported models or you can get erroneous results.

Value

None

Author(s)

Jeff Laake

See Also

import.chdata

Examples

data(dipper)
dipper$numeric.sex=as.numeric(dipper$sex)-1
dipper.processed=process.data(dipper,group="sex")
export.chdata(dipper.processed, filename="dipper", 
          covariates="numeric.sex",replace=TRUE)
file.remove("dipper.inp")

#
# Had sex been used in place of numeric.sex in the above command, 
# MARK would have been unable to use it as a covariate
# because it is not a numeric field

Description

Creates a .Rinp, .inp and optionally renamed output files that can be imported into MARK to create a new MARK project with all the data and output files.

Usage

export.MARK(
  x,
  project.name,
  model = NULL,
  replace = FALSE,
  chat = 1,
  title = "",
  ind.covariates = "all"
)

Arguments

Details

If you use Nest model and NestAge covariate you should use the processed data list (model$data) from a model using NestAge in the formula because the necessary individual covariates are added to the processed data list. Also, use default of ind.covariates="all".

After running this function to export the data and models to the working directory, start program MARK and select the File/RMARK Import menu item. Navigate to the working directory and select the "project.name".Rinp file. MARK will take over and create the project files and will import any specified model output files.

DO NOT use a "project name" that is the same as a filename you use as input (*.inp) to RMark because this function will create a file "project name".inp and it would over-write your existing file and the format could change. If you try this, the code will return an error that the filename is invalid because the file already exists.

Value

None

Author(s)

Jeff Laake

See Also

export.chdata, export.model

Examples

# This example is excluded from testing to reduce package check time
data(mallard)
Dot=mark(mallard,nocc=90,model="Nest",
	model.parameters=list(S=list(formula=~1)),delete=TRUE)
mallard.proc=process.data(mallard,nocc=90,model="Nest")
# removed # to use export.MARK
#export.MARK(mallard.proc,"mallard",Dot,replace=TRUE)
data(robust)
time.intervals=c(0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0)
S.time=list(formula=~time)
p.time.session=list(formula=~-1+session:time,share=TRUE)
GammaDoublePrime.random=list(formula=~time,share=TRUE)
model.1=mark(data = robust, model = "Robust",
	time.intervals=time.intervals,
	model.parameters=list(S=S.time,
			GammaDoublePrime=GammaDoublePrime.random,p=p.time.session),delete=TRUE)
robust.proc=process.data(data = robust, model = "Robust",
	time.intervals=time.intervals)
#remove # to use export.MARK
#export.MARK(robust.proc,"robust",	model.1,replace=TRUE)

Description

Creates renamed versions of the output,vcv and residual files so they can be appended into a MARK .dbf file.

Usage

export.model(model, replace = FALSE)

Arguments

Details

If model is a marklist then it exports each model in the marklist. The function simply copies the files with new names so the MARK interface will recognize them. The marknnn.out is copied as marknnnY.tmp, marknnn.res is copied as marknnnx.tmp and marknnn.vcv is copied as marknnnV.tmp. You can create a MARK .dbf by using export.chdata to create an input file for MARK, opening MARK (MARKINT.EXE) to create a new .dbf with the input file, and then using the Output/Append to select the output file (marknnnY.tmp) to append the model with its files. Then you can use any facilities of MARK that are not already included in RMark.

***Warning*** Make sure that you use the .inp created by export.chdata with your processed data to create the MARK .dbf file rather than using a separate similar .inp file. It is essential that the group structure and ordering of groups matches between the .inp file and the exported models or you can get erroneous results.

Value

None

Author(s)

Jeff Laake

See Also

export.chdata

Examples

data(dipper)
mymodel=mark(dipper,threads=1,delete=TRUE)
# remove # to use export.model
#export.model(mymodel,replace=TRUE)

Description

Miscellaneous set of functions that can be used with results from the package.

Usage

extract.indices(model,parameter,df)
 
        nat.surv(model,df)

         pop.est(ns,ps,design,p.vcv)

         compute.Sn(x,df,criterion)

         logitCI(x,se)

         search.output.files(x,string)

Arguments

Details

Function extract.indices extracts the parameter indices from the parameter index matrices (PIMS) for a particular type of parameter that match a set of group numbers and rows and columns that are defined in the dataframe df. It returns a vector of indices which can be used to specify the set of real parameters to be extracted by covariate.predictions using the index column in data or the indices argument. If df is NULL, it returns a dataframe with all of the indices with model.index being the unique index across all parameters and the par.index which is an index to the row in the design data. If parameter is NULL then the the dataframe is given for all of the parameters.

Function nat.surv produces estimates of natural survival (Sn) from total survival (S) and recovery rate (r) from a joint live-dead model in which all harvest recoveries are reported. In that case, Taylor et al 2005 suggest the following estimator of natural survival Sn=S + (1-S)*r. The arguments for the function are a mark model object and a dataframe df that defines the set of groups and times (row,col) for the natural survival computations. It returns a list with elements: 1) Sn - a vector of estimates for natural survival; one for each entry in df and 2) vcv - a variance-covariance matrix for the estimates of natural survival.

Function pop.est produces estimates of abundance using a vector of counts of animals captured (ns) and estimates of capture probabilities (ps). The estimates can be aggregated or averaged using the design matrix argument. If individual estimates are needed, use an nxn identity matrix for design where n is the length of ns. To get a total of all the estimates use a nx1 column matrix of 1s. Any other design matrix can be specified to subset, aggregate and/or average the estimates. The argument p.vcv is needed to compute the variance-covariance matrix for the abundance estimates using the formula described in Taylor et al. (2002). The function returns a list with elements: 1) Nhat - a vector of abundance estimates and 2) vcv - variance-covariance matrix for the abundance estimates.

Function Compute.Sn creates list structure for natural survival using nat.surv to be used for model averaging natural survival estimates (e.g., model.average(compute.Sn(x,df,criterion))). It returns a list with elements estimates, vcv, weight: 1) estimates - matrix of estimates of natural survival, 2)vcv - list of var-cov matrix for the estimates, and 3) weight - vector of model weights.

Function search.output.filessearches for occurrence of a specific string in output files associated with models in a marklist x. It returns a vector of model numbers in the marklist which have an output file containing the string.

Author(s)

Jeff Laake

References

TAYLOR, M. K., J. LAAKE, H. D. CLUFF, M. RAMSAY and F. MESSIER. 2002. Managing the risk from hunting for the Viscount Melville Sound polar bear population. Ursus 13: 185-202.

TAYLOR, M. K., J. LAAKE, P. D. MCLOUGHLIN, E. W. BORN, H. D. CLUFF, S. H. FERGUSON, A. ROSING-ASVID, R. SCHWEINSBURG and F. MESSIER. 2005. Demography and viability of a hunted population of polar bears. Arctic 58: 203-214.

Examples

# This example is excluded from testing to reduce package check time
# Example of computing N-hat for occasions 2 to 7 for the p=~time model
data(dipper)
md=mark(dipper,model.parameters=list(p=list(formula=~time),
       Phi=list(formula=~1)),delete=TRUE)
# Create a matrix from the capture history strings 
xmat=matrix(as.numeric(unlist(strsplit(dipper$ch,""))),
      ncol=nchar(dipper$ch[1]),byrow=TRUE)
# sum number of captures in each column but don't use the first 
# column because p[1] can't be estimated
ns=colSums(xmat)[-1]
# extract the indices and then get covariate predictions for p(2),...,p(7)
# which are row-colums 1-6 in PIM for p 
p.indices=extract.indices(md,"p",df=data.frame(group=rep(1,6),
   row=1:6,col=1:6))
p.list=covariate.predictions(md,data=data.frame(index=p.indices))
# call pop.est using diagonal design matrix to get 
# separate estimate for each occasion
pop.est(ns,p.list$estimates$estimate,
  design=diag(1,ncol=6,nrow=6),p.list$vcv)

Description

Extracts the lnl, AICc, npar, beta and real estimates and returns a list of these results for inclusion in the mark object. The elements beta and real are dataframes with fields estimate,se,lcl,ucl. This function was written for internal use and is called by run.mark.model. It is documented here for more advanced users that might want to modify the code or adapt for their own use.

Usage

extract.mark.output(out, model, adjust, realvcv = FALSE, vcvfile)

Arguments

Value

result: list of extracted output elements

Author(s)

Jeff Laake

See Also

run.mark.model

Description

Replaces covariate names in design matrix with specific values to compute estimates of real parameters at those values using the dataframe from find.covariates after any value replacement.

Usage

fill.covariates(model, values)

Arguments

Details

The design matrix for a MARK model with individual covariates contains the covariate names used in the model. In computing the real parameters for the encounter history of an individual it replaces instances of covariate names with the individual covariate values. This function replaces the cells in the design matrix that contain individidual covariates with user-specified values which is an edited version (if needed) of the dataframe returned by find.covariates.

Value

New design matrix with user-defined covariate values entered in place of covariate names

Author(s)

Jeff Laake

See Also

find.covariates, compute.real

Examples

data(dipper)
dipper$nsex=as.numeric(dipper$sex)-1
dipper$weight=rnorm(294)   
#NOTE:  This generates random valules for the weights so the answers using 
# ~weight will vary each time it is run
mod=mark(dipper,model.parameters=list(Phi=list(formula=~nsex+weight)),delete=TRUE)
# Show approach using individual calls to find.covariates, fill.covariates 
# and compute.real
fc=find.covariates(mod,dipper)
fc$value[fc$var=="nsex"]=0 # assign sex value to Female
design=fill.covariates(mod,fc) # fill design matrix with values
# compute and output survivals for females at average weight
female.survival=compute.real(mod,design=design)[1,] 
female.survival
# Next show same thing with a call to compute.real and a data frame for 
# females and then males
# compute and output survivals for females at average weight
female.survival=compute.real(mod,data=
      data.frame(nsex=0,weight=mean(dipper$weight)))[1,] 
female.survival
male.survival=compute.real(mod,data=data.frame(nsex=1,
         weight=mean(dipper$weight)))[1,] 
male.survival
# Fit model using sex as a group/factor variable and 
# compute v-c matrix for estimates
mod=mark(dipper,groups="sex",
     model.parameters=list(Phi=list(formula=~sex+weight)),delete=TRUE)
survival.by.sex=compute.real(mod,data=dipper,vcv=TRUE)
survival.by.sex$real[1:2]  # estimates
survival.by.sex$se.real[1:2] # std errors
survival.by.sex$vcv.real[1:2,1:2] # v-c matrix
survival.by.sex$vcv.real[1,2]/prod(survival.by.sex$se.real[1:2]) 
# sampling correlation of the estimates