Lab 02 - Analysis in RMark
WILD 502 - Jay Rotella
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.ziplibrary(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.3074218Look 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.6540015Work 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 maleWorking 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)