--- title: "Lab 08 - Analyses in RMark" author: "WILD 502 - Jay Rotella" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` This markdown document runs the POPAN models for data on Coho salmon that you analyzed directly in MARK in lab 8. Here, we're dealing with 4 types of parameters: 'Phi', 'p', 'pent', and 'N'. ### Bring in the Data ```{r, warning=FALSE} library(RMark) coho <- convert.inp("http://www.montana.edu/rotella/documents/502/chaseAD.inp", group.df = NULL) head(coho) ``` ### Process the data ```{r} # Process data # NOTE: for this data set, we need to provide the time # intervals because they are not all the same length coho.pr <- process.data(coho, begin.time = 1, model = "POPAN", time.intervals = c(1.5, 1, 1, 1, 1, 1, 1.5)) # # Create default design data coho.ddl = make.design.data(coho.pr) # examine the design data coho.ddl ``` ### Build function for creating models Here, we set up a function that contains the structures for 'Phi', 'p', 'pent' and 'N'. It will run all combinations of the structures for each parameter type when executed. ```{r} run.coho <- function() { # # Define range of models for each parameter type. Here, # we only have a few models as the labs focused more on how # how work with the output than on setting up a wide variety # of competing models # Define model for Phi Phi.time = list(formula = ~ time) # Define range of models for p p.dot = list(formula = ~ 1) p.time = list(formula = ~ time) # Define model for pent (probability of entry) pent.time = list(formula = ~ time) # Define model for N N.dot = list(formula = ~ 1) # Run assortment of models Phi.time.p.dot.pent.time.N = mark( coho.pr, coho.ddl, model.parameters = list( Phi = Phi.time, p = p.dot, pent = pent.time, N = N.dot)) # Note this initial run will penalize for too many parameters because # we are not adjusting the parameter count for confounded parameters Phi.time.p.time.pent.time.N = mark( coho.pr, coho.ddl, model.parameters = list( Phi = Phi.time, p = p.time, pent = pent.time, N = N.dot)) # Here, we adjust the parameter count to the correct value Phi.time.p.time.pent.time.N = adjust.parameter.count(Phi.time.p.time.pent.time.N, 21) # Return model table and list of models return(collect.models()) } ``` ### Run the Models and Examine the Output It's very simple to then run the function and each of the models. As implemented here, the output from each of the models appears in the console where it can be reviewed when the function written above is called. If you don't want the output from every model to appear automatically, see the note above the "mark.wrapper" command for how to set the option of suppressing the output. One can also examine model-specific output in other ways as shown below. ```{r} coho.results <- run.coho() ``` ### Examine Model-Selection Table Once the models are run, we can examine the model-selection table. We can also examine model-specific output. ```{r} options(width = 150) coho.results ``` ### Look at output from top model ```{r} # View the model names names(coho.results) # examine the output from top-ranked model (#1) and # store top model in object with short name top <- coho.results$Phi.time.p.dot.pent.time.N # Check which link functions were used for the parameters # Notice that RMark used parameter-specific link functions that # make sense here, i.e., the logit link for Phi & p, the mlogit for pent, # and a log link for N. top$links # look at summary of top model's output summary(top, showall = FALSE) # store and examine estimates of Phi, p, pent, and N # First examine estimates of the real parameters # to see how results are stored. We drop 'fixed' & 'note' columns, so # we can use the 'round' command round(top$results$real[ , 1:4], 5) # for Phi in top model there are 7 estimates top.Phi <- top$results$real[1:7, ] top.Phi # for p in top model there is 1 estimate # and it's in the 8th row of output top.p <- top$results$real[8, ] top.p # Store and examine the estimates of pent top.pent <- top$results$real[10:16, ] round(top.pent[ ,1:4], 5) # Store estimate for N top.N <- top$results$real[9, ] top.N ``` ### Calculate the estimated number of births To calculate the estimates of B[i], you need to multiply N by b[i]. Remember, to get b[0], you subtract the sum of all the b[i] for i = 1, 2, ..., 7 from 1. ```{r} # calculate the proportion of the population present at the start. # To do this we need the b[i] for each occasion b <- top.pent$estimate b0 <- 1 - sum(top.pent$estimate) # # Examine the b[i] including b[0] round(c(b0, b), 4) # Next, calculate the net births by occasion B.occ <- b*top.N$estimate round(B.occ, 4) # B[0] is simply the proportion of the superpopulation that was # present at the start of the study (b[0]*N-hat) (B0 <- b0*top.N$estimate) # Sum up all the B[i] including B0 sum(B0, B.occ) ``` ### Work out N[i] values For this, we have to deal with the losses on capture and the fact that not all time intervals are the same between the various occasions. The details are provided on page 12-20 of the C&W chapter. ```{r} # review values for B[i] round(c(B0, B.occ), 4) # calculate and store values of N[i] for each occasion iteratively # set up vector for containing N by occ N.occ <- vector("numeric", 8) # Store starting number N.occ[1] <- B0 # Store time intervals for use in calculations time.intervals <- c(1.5, 1, 1, 1, 1, 1, 1.5) # Store the losses on capture (see Table 12.6 on page 12-15 of C&W) losses <- c(0, 1, 11, 2, 8, 8, 6, 10) # Calculate the N.occ values for (i in 2:8) { N.occ[i] = (N.occ[i - 1] - losses[i - 1]) * (top.Phi$estimate[i - 1]^time.intervals[i - 1]) + B.occ[i - 1] } # examine estimated numbers round(cbind(N.occ = N.occ, B.occ = c(B.occ, NA)), 4) ``` ### Further work that can be done Chapter 12 goes through more modeling that can be done including additional models of the adults and combined modeling for adults and jacks. Also, you could work with the delta method to obtain estimates of the SE's for the quantities we derived above for B.occ and N.occ.