One of the most frequent questions I get is "why do the results from MARK and RMark differ"? Initially it was so frequent that Evan created the locked posting "when RMark and MARK results (seem to) differ | read first".
The content of that posting is correct but I want to revisit and elaborate some. First of all, I'll never discourage anyone from comparing the results of MARK and RMark, just remember that RMark doesn't fit models, MARK does and if there are differences it is for one of the reasons listed in that posting. In fact, I encourage you to compare the results because you'll become more comfortable with RMark and it helps me find any bugs that many still be lurking in RMark. My recent exchange with Kristin helped me discover a bug in model.average.list in which the non-revised model average variance formula was being used if the full v-c matrix was not specified. She was using it to model average derived parameters and noted the difference.
The other differences she noted were quite small ones (1e-5 to 1e-8) and resulted from differences in the design matrix formulation. Although in theory it shouldn't make a difference, you may get slight differences when the design matrices differ because the model fitting algorithm is numerical. This has been noted several times on the list and it primarily occurs with small data sets. One of the primary DM differences is in regard to factor variables like time. RMark uses the R contrast setting which by default is the treatment contrast which uses the first factor level as the intercept and in MARK Gary uses the SAS convention of the last factor level being the intercept. In the example below, I show how you can modify the level that is used as the intercept with the relevel function on the design data. Also, in the example I use I(1-c) instead of c in the formula to match the DM Kristin created in MARK. BOTTOM LINE: same model can be specified with many different DMs and these can result in very small differences in the values but if you construct the same DM with a formula that you did in the MARK interface you WILL get the same answer because it is the same code fitting the model. The best way to assure yourself of this is to use export.MARK to export the models you fitted into MARK and then re-run them in the MARK interface.
The deviance value put in the marklist model table is the residual deviance. If you are comparing models with and without covariates the deviance values will be quite different because the null deviance with
covariates is set to 0 by MARK. You can use a switch in model.table to use -2lnl instead of deviance.
regards--jeff
- Code: Select all
#convert INP file
sq=convert.inp("AG2011AUG.inp", use.comments=TRUE)
#step one: process.data
sq.HFH=process.data(sq,model="HugFullHet")
sq.Hug=process.data(sq,model="Huggins")
#step two: make.design.data
sq.HFHddl=make.design.data(sq.HFH)
sq.Hugddl=make.design.data(sq.Hug)
# change level for time and mixture factor variable
sq.HFHddl$p$time=relevel(sq.HFHddl$p$time,"6")
sq.HFHddl$c$time=relevel(sq.HFHddl$c$time,"6")
sq.Hugddl$p$time=relevel(sq.Hugddl$p$time,"6")
sq.Hugddl$c$time=relevel(sq.Hugddl$c$time,"6")
sq.HFHddl$p$mixture=relevel(sq.HFHddl$p$mixture,"2")
sq.HFHddl$c$mixture=relevel(sq.HFHddl$c$mixture,"2")
### MODEL SPECIFICATIONS ####3
# p(~1)c() or M0
sq.HugFullHet.M0=mark(sq.Hug,sq.Hugddl,model.parameters=list(p=list(formula=~1,share=TRUE)))
sq.HugFullHet.M0.alt=mark(sq.HFH,sq.HFHddl,model.parameters=list(p=list(formula=~1,share=TRUE),pi=list(formula=~1,fixed=1))
# p(t)c() or Mt
sq.HugFullHet.Mt=mark(sq.Hug,sq.Hugddl,model.parameters=list(p=list(formula=~time,share=TRUE)))
# p(~1)c(~1) or Mb
sq.HugFullHet.Mb=mark(sq.Hug,sq.Hugddl,model.parameters=list(p=list(formula=~1,share=FALSE),c=list(formula=~1)))
# p(~1)c(),mixture or Mh
sq.HugFullHet.Mh=mark(sq.HFH,sq.HFHddl,model.parameters=list(p=list(formula=~mixture,share=TRUE)))
# p(t)c(t) or Mtb
sq.HugFullHet.Mtb=mark(sq.Hug,sq.Hugddl,model.parameters=list(p=list(formula=~time+I(1-c),share=TRUE)))
# p(~1)c(~1),mixture or Mbh
sq.HugFullHet.Mbh=mark(sq.HFH,sq.HFHddl,,model.parameters=list(p=list(formula=~mixture+I(1-c),share=TRUE)))
# p(t)c(),mixture or Mth
sq.HugFullHet.Mth=mark(sq.HFH,sq.HFHddl,,model.parameters=list(p=list(formula=~time+mixture,share=TRUE)))
# p(t)c(t),mixture or Mtbh
sq.HugFullHet.Mtbh=mark(sq.HFH,sq.HFHddl,,model.parameters=list(p=list(formula=~time+mixture+I(1-c),share=TRUE)))
# To check out the design matrix to see if it looks like the same design as in MARK interface
#sq.HugFullHet.M0$design.matrix
#sq.HugFullHet.Mb$design.matrix
#sq.HugFullHet.Mt$design.matrix
#sq.HugFullHet.Mh$design.matrix
#sq.HugFullHet.Mtb$design.matrix
#sq.HugFullHet.Mth$design.matrix
#sq.HugFullHet.Mbh$design.matrix
#sq.HugFullHet.Mtbh$design.matrix
### SEE DERIVED N-hat RESULTS #####3
# If you want to see more digits use; I think Gary's output is truncating and R is rounding
options(digits=8)
sq.HugFullHet.M0$results$derived
sq.HugFullHet.Mb$results$derived
sq.HugFullHet.Mt$results$derived
sq.HugFullHet.Mh$results$derived
sq.HugFullHet.Mth$results$derived
sq.HugFullHet.Mtb$results$derived
sq.HugFullHet.Mbh$results$derived
sq.HugFullHet.Mtbh$results$derived