### 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)