Introduction

In lab this week, we have known fate data for 1 interval of time and the input file also contains 4 individual covariates (Area, Sex, Mass, & Length). Here, we’ll use RMark and work with the same models that you evaluated directly in MARK for lab 2.

Import and process the data

library(RMark)
## This is RMark 2.2.4
##  Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
library(ggplot2)
# bring in the data file
fawns <- convert.inp("http://www.montana.edu/rotella/documents/502/lab02-fawns.inp",
                     covariates = c("area", "sex", "mass", "length"))

# add "length.sq" variable for use in quadratic models
fawns$length.sq <- (fawns$length)^2
head(fawns)
##   ch freq area sex mass length length.sq
## 1 11    1    1   0 34.0  126.0  15876.00
## 2 11    1    1   0 33.4  128.5  16512.25
## 3 11    1    1   0 29.9  130.0  16900.00
## 4 11    1    1   0 31.0  113.0  12769.00
## 5 11    1    1   0 26.7  116.0  13456.00
## 6 11    1    1   0 23.9  124.0  15376.00
# convert area & sex variables to factors so that we can treat them as 
# grouping variables in the analysis to ease working with output later
fawns$area <- factor(fawns$area, 
                     levels = c(0, 1), 
                     labels = c("control", "treatment"))
fawns$sex  <- factor(fawns$sex,
                     levels = c(0, 1),
                     labels = c("female", "male"))
# Store versions of fate for use in final plotting
# first a factor version that has nice labels
fawns$fate <- factor(as.numeric(fawns$ch), 
                     levels = c(11, 10),
                     labels = c("died", "lived"))
# second, a numeric version that's 0 = dead, 1 = lived
fawns$Fate <- 2 - as.numeric(fawns$fate)

summary(fawns)
##       ch                 freq          area        sex          mass      
##  Length:115         Min.   :1   control  :59   female:57   Min.   :22.80  
##  Class :character   1st Qu.:1   treatment:56   male  :58   1st Qu.:31.90  
##  Mode  :character   Median :1                              Median :33.60  
##                     Mean   :1                              Mean   :34.38  
##                     3rd Qu.:1                              3rd Qu.:36.60  
##                     Max.   :1                              Max.   :72.00  
##      length        length.sq        fate         Fate       
##  Min.   :108.0   Min.   :11664   died :68   Min.   :0.0000  
##  1st Qu.:120.0   1st Qu.:14400   lived:47   1st Qu.:0.0000  
##  Median :124.0   Median :15376              Median :1.0000  
##  Mean   :123.2   Mean   :15217              Mean   :0.5913  
##  3rd Qu.:127.0   3rd Qu.:16129              3rd Qu.:1.0000  
##  Max.   :135.5   Max.   :18360              Max.   :1.0000
# we know from other work on this dataset that there is 1 record with 
#  a very high value for mass, which we want to drop
fawns2 = fawns[fawns$mass < 70, ]

Process the data for use in RMark

# Process data
fawns.processed = process.data(fawns2, 
                               model = "Known", 
                               groups = c("area", "sex"))
#
# Create default design data
fawns.ddl = make.design.data(fawns.processed)

Build function for creating models

Here, we write a function that establishes which models are to be run in MARK.

# setup a function 
run.fawns <- function() {
#  Define range of models for S
  S.dot = list(formula =  ~ 1)
  S.area = list(formula =  ~ area)
  S.mass = list(formula =  ~ mass)
  S.length = list(formula =  ~ length)
  S.length.sq = list(formula =  ~ length + length.sq)
  S.sex = list(formula =  ~ sex)
  S.area.length = list(formula =  ~ area + length)
  S.length.mass = list(formula =  ~ length + mass)
  # sex-specific intercepts & common slope
  S.sex.length = list(formula =  ~ sex + length)
  # 1 intercept & sex-specific slopes
  S.length.by.sex = list(formula =  ~ sex:length)
  # sex-specific intercepts & slopes
  S.sex.x.length = list(formula =  ~ sex * length)

# Create model list
  model.list = create.model.list("Known")

  # NOTE: to avoid having all the output for each model appear when you
  # call the function, add ', output=FALSE' after 'ddl=fawns.ddl' below.
  # Here, I don't do that so you can see the output for each model,
  # but this might not be desired if you have many models.
  fawns.results = mark.wrapper(model.list, data = fawns.processed,
    ddl = fawns.ddl, invisible = TRUE, threads = 2)
    
  # Return model table and list of models
  return(fawns.results)
}

Run the models

It’s worth noting that here, the output from each of the models appears in the console where it can be reviewed. One can also examine model-specific output in other ways as shown below. And, if you don’t want to see the output from every model appear in the output here, see the note above for the use of the “mark.wrapper” command where one has the option of suppressing the output.

fawns.results = run.fawns()
## 
## S.area
## 
## Output summary for Known model
## Name : S(~area) 
## 
## Npar :  2
## -2lnL:  154.39
## AICc :  158.4981
## 
## Beta
##                   estimate        se        lcl       ucl
## S:(Intercept)   -0.4198538 0.2684207 -0.9459585 0.1062508
## S:areatreatment  0.1321718 0.3807444 -0.6140873 0.8784309
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.3965517
## Group:areatreatment.sexfemale 0.4285714
## Group:areacontrol.sexmale     0.3965517
## Group:areatreatment.sexmale   0.4285714
## 
## S.area.length
## 
## Output summary for Known model
## Name : S(~area + length) 
## 
## Npar :  3
## -2lnL:  149.9375
## AICc :  156.1557
## 
## Beta
##                    estimate        se         lcl        ucl
## S:(Intercept)   -10.1984070 4.7980557 -19.6025960 -0.7942173
## S:areatreatment   0.0712668 0.3894067  -0.6919703  0.8345039
## S:length          0.0794145 0.0388388   0.0032904  0.1555386
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.3998622
## Group:areatreatment.sexfemale 0.4170799
## Group:areacontrol.sexmale     0.3998622
## Group:areatreatment.sexmale   0.4170799
## 
## S.dot
## 
## Output summary for Known model
## Name : S(~1) 
## 
## Npar :  1
## -2lnL:  154.5106
## AICc :  156.5463
## 
## Beta
##                estimate        se        lcl       ucl
## S:(Intercept) -0.354545 0.1902681 -0.7274706 0.0183806
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.4122807
## Group:areatreatment.sexfemale 0.4122807
## Group:areacontrol.sexmale     0.4122807
## Group:areatreatment.sexmale   0.4122807
## 
## S.length
## 
## Output summary for Known model
## Name : S(~length) 
## 
## Npar :  2
## -2lnL:  149.971
## AICc :  154.0791
## 
## Beta
##                  estimate        se         lcl        ucl
## S:(Intercept) -10.2340400 4.7978013 -19.6377300 -0.8303488
## S:length        0.0799874 0.0387469   0.0040435  0.1559314
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.4082924
## Group:areatreatment.sexfemale 0.4082924
## Group:areacontrol.sexmale     0.4082924
## Group:areatreatment.sexmale   0.4082924
## 
## S.length.by.sex
## 
## Output summary for Known model
## Name : S(~sex:length) 
## 
## Npar :  3
## -2lnL:  142.8986
## AICc :  149.1168
## 
## Beta
##                       estimate        se         lcl        ucl
## S:(Intercept)      -13.2950770 5.1909361 -23.4693120 -3.1208418
## S:sexfemale:length   0.1090090 0.0423959   0.0259130  0.1921050
## S:sexmale:length     0.1002673 0.0414444   0.0190362  0.1814984
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.5365586
## Group:areatreatment.sexfemale 0.5365586
## Group:areacontrol.sexmale     0.2826384
## Group:areatreatment.sexmale   0.2826384
## 
## S.length.mass
## 
## Output summary for Known model
## Name : S(~length + mass) 
## 
## Npar :  3
## -2lnL:  149.9709
## AICc :  156.189
## 
## Beta
##                    estimate        se         lcl       ucl
## S:(Intercept) -1.026642e+01 5.3917051 -20.8341610 0.3013239
## S:length       8.048000e-02 0.0538664  -0.0250980 0.1860581
## S:mass        -8.335333e-04 0.0633061  -0.1249135 0.1232464
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.4082881
## Group:areatreatment.sexfemale 0.4082881
## Group:areacontrol.sexmale     0.4082881
## Group:areatreatment.sexmale   0.4082881
## 
## S.length.sq
## 
## Output summary for Known model
## Name : S(~length + length.sq) 
## 
## Npar :  3
## -2lnL:  145.888
## AICc :  152.1062
## 
## Beta
##                   estimate           se          lcl          ucl
## S:(Intercept) -212.9691700 6.3699964000 -225.4543600 -200.4839800
## S:length         3.3680493 0.0932691000    3.1852419    3.5508567
## S:length.sq     -0.0133118 0.0004137664   -0.0141228   -0.0125008
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.3933261
## Group:areatreatment.sexfemale 0.3933261
## Group:areacontrol.sexmale     0.3933261
## Group:areatreatment.sexmale   0.3933261
## 
## S.mass
## 
## Output summary for Known model
## Name : S(~mass) 
## 
## Npar :  2
## -2lnL:  152.281
## AICc :  156.3891
## 
## Beta
##                 estimate        se        lcl       ucl
## S:(Intercept) -2.6154908 1.5533130 -5.6599843 0.4290027
## S:mass         0.0661987 0.0450065 -0.0220142 0.1544115
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.4106284
## Group:areatreatment.sexfemale 0.4106284
## Group:areacontrol.sexmale     0.4106284
## Group:areatreatment.sexmale   0.4106284
## 
## S.sex
## 
## Output summary for Known model
## Name : S(~sex) 
## 
## Npar :  2
## -2lnL:  150.0979
## AICc :  154.206
## 
## Beta
##                 estimate        se        lcl        ucl
## S:(Intercept)  0.0350912 0.2649472 -0.4842054  0.5543878
## S:sexmale     -0.8082811 0.3890933 -1.5709040 -0.0456582
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.5087719
## Group:areatreatment.sexfemale 0.5087719
## Group:areacontrol.sexmale     0.3157895
## Group:areatreatment.sexmale   0.3157895
## 
## S.sex.length
## 
## Output summary for Known model
## Name : S(~sex + length) 
## 
## Npar :  3
## -2lnL:  143.2472
## AICc :  149.4654
## 
## Beta
##                  estimate        se         lcl        ucl
## S:(Intercept) -12.6958420 5.1232179 -22.7373490 -2.6543346
## S:sexmale      -1.0554989 0.4181604  -1.8750932 -0.2359046
## S:length        0.1040179 0.0417841   0.0221212  0.1859147
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.5325286
## Group:areatreatment.sexfemale 0.5325286
## Group:areacontrol.sexmale     0.2838994
## Group:areatreatment.sexmale   0.2838994
## 
## S.sex.x.length
## 
## Output summary for Known model
## Name : S(~sex * length) 
## 
## Npar :  4
## -2lnL:  140.1658
## AICc :  148.5328
## 
## Beta
##                     estimate         se         lcl        ucl
## S:(Intercept)    -22.7241270  8.4186677 -39.2247160 -6.2235377
## S:sexmale         17.5796290 10.8831710  -3.7513878 38.9106450
## S:length           0.1858880  0.0686609   0.0513127  0.3204633
## S:sexmale:length  -0.1507537  0.0881650  -0.3235571  0.0220496
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.5491318
## Group:areatreatment.sexfemale 0.5491318
## Group:areacontrol.sexmale     0.3074218
## Group:areatreatment.sexmale   0.3074218

Look at the model-selection table

fawns.results
##                     model npar     AICc DeltaAICc      weight    Deviance
## 11       S(~sex * length)    4 148.5328 0.0000000 0.363470826 140.1657900
## 5          S(~sex:length)    3 149.1168 0.5840393 0.271423238 142.8986200
## 10       S(~sex + length)    3 149.4654 0.9326393 0.228007692 143.2472200
## 7  S(~length + length.sq)    3 152.1062 3.5733993 0.060885763 145.8879800
## 4              S(~length)    2 154.0791 5.5463756 0.022703321 149.9710300
## 9                 S(~sex)    2 154.2060 5.6732356 0.021307971   0.3618895
## 2       S(~area + length)    3 156.1557 7.6229593 0.008038314 149.9375400
## 6       S(~length + mass)    3 156.1890 7.6562693 0.007905545 149.9708500
## 8                S(~mass)    2 156.3891 7.8563056 0.007153103 152.2809600
## 3                   S(~1)    1 156.5463 8.0135118 0.006612377   4.7745635
## 1                S(~area)    2 158.4981 9.9653456 0.002491852   4.6540015

Work with the top model

You would not need to re-run the top model since we’ve already run it. But, by doing it this way, you can see how you can to run a model of interest without using a function as was done earlier.

# here's how to run a model outside the function
top <- mark(fawns.processed,ddl = fawns.ddl, 
            model = "Known",
            model.parameters = list("S" = list(formula =  ~ sex * length)))
## 
## Output summary for Known model
## Name : S(~sex * length) 
## 
## Npar :  4
## -2lnL:  140.1658
## AICc :  148.5328
## 
## Beta
##                     estimate         se         lcl        ucl
## S:(Intercept)    -22.7241260  8.4186686 -39.2247170 -6.2235355
## S:sexmale         17.5796310 10.8831740  -3.7513892 38.9106520
## S:length           0.1858880  0.0686609   0.0513126  0.3204633
## S:sexmale:length  -0.1507537  0.0881650  -0.3235571  0.0220496
## 
## 
## Real Parameter S
##                                       2
## Group:areacontrol.sexfemale   0.5491318
## Group:areatreatment.sexfemale 0.5491318
## Group:areacontrol.sexmale     0.3074218
## Group:areatreatment.sexmale   0.3074218
# look at beta-hats
top$results$beta
##                     estimate         se         lcl        ucl
## S:(Intercept)    -22.7241260  8.4186686 -39.2247170 -6.2235355
## S:sexmale         17.5796310 10.8831740  -3.7513892 38.9106520
## S:length           0.1858880  0.0686609   0.0513126  0.3204633
## S:sexmale:length  -0.1507537  0.0881650  -0.3235571  0.0220496
# create values of length to use for predictions
min.length = min(fawns2$length)
max.length = max(fawns2$length)
length.values = seq(from = min.length, to = max.length, length = 100)

# determine which parameter indices go with males and females
fawns.ddl
## $S
##   par.index model.index           group age time Age Time      area    sex
## 1         1           1   controlfemale   1    2   1    0   control female
## 2         2           2 treatmentfemale   1    2   1    0 treatment female
## 3         3           3     controlmale   1    2   1    0   control   male
## 4         4           4   treatmentmale   1    2   1    0 treatment   male
## 
## $pimtypes
## $pimtypes$S
## $pimtypes$S$pim.type
## [1] "all"
# make predictions for rows of 'fawns.ddl' associated with
# females (par.index=1) and males (par.index=3)
pred.top <- covariate.predictions(top, 
              data = data.frame(length = length.values, 
                                length.sq = length.values ^ 2),
              indices = c(1, 3))

# store values of sex in pred.top
female.rows <- which(pred.top$estimates$par.index == 1)
male.rows <- which(pred.top$estimates$par.index == 3)
pred.top$estimates$sex <- NA
pred.top$estimates$sex[female.rows] <- "female"
pred.top$estimates$sex[male.rows] <- "male"
head(pred.top$estimates)
##   vcv.index model.index par.index   length length.sq   estimate         se
## 1         1           1         1 108.0000  11664.00 0.06609868 0.06414453
## 2         2           2         3 108.0000  11664.00 0.20587191 0.15729338
## 3         3           1         1 108.2778  11724.08 0.06935843 0.06589106
## 4         4           2         3 108.2778  11724.08 0.20747206 0.15578488
## 5         5           1         1 108.5556  11784.31 0.07276642 0.06764200
## 6         6           2         3 108.5556  11784.31 0.20908138 0.15425279
##           lcl       ucl fixed    sex
## 1 0.009149198 0.3517066       female
## 2 0.037843503 0.6308201         male
## 3 0.009977715 0.3553053       female
## 4 0.039270820 0.6263872         male
## 5 0.010879836 0.3589333       female
## 6 0.040747949 0.6219434         male

Working with output: Plotting

Here, we plot estimates of “S” for male and female fawns of different lengths. We also want to plot CI’s for the values. To do so, we use the predicted values we created in the last step. Here, I use the ggplot2 package, which makes nice graphs.

pred <- pred.top$estimates
# build and store the plot in object 'p'
p <- ggplot(pred, aes(x = length, y = estimate, group = sex)) +
  geom_line(size = 1.5, aes(colour = sex)) +
  geom_ribbon(aes(ymin = lcl, ymax = ucl), alpha = 0.2) +
  scale_colour_brewer(palette = "Set1") +
  theme(legend.position = c(0.15, 0.8),
  legend.justification = c(1, 0)) +
    xlab("Body Length (cm)") + 
    ylab("Estimated Survival Rate") + 
  geom_jitter(data = fawns2, aes(y = Fate, color = sex),
              height = 0.025)
# print the plot
p

# if you're running this code, it's wise to cleanup files
# you can use the next 2 lines without comments to do that
#rm(list=ls(all=TRUE))
#cleanup(ask = FALSE)