splines in RMark

[jlaake]

I wanted to pass on a modelling approach that was suggested to me by Alexey Altukhov via Eli Gurarie to fit splines with RMark. Below is an example with the dipper data

Code: Select all

library(RMark)
library(splines)
library(ggplot2)
data(dipper)
dipper.proc=process.data(dipper[dipper$sex=="Female",],model="CJS")
dipper.ddl=make.design.data(dipper.proc)
mod.Phis.pt=mark(dipper.proc,dipper.ddl,model.parameters=list(Phi=list(formula=~bs(Time))))
predPhi <- summary(mod.Phis.pt, se=TRUE)$real$Phi[1:6,]
pred.sm <- data.frame(spline(predPhi$Time, predPhi$estimate, n=100))
pred.sm$lcl <- spline(predPhi$Time, predPhi$lcl, n=100)$y
pred.sm$ucl <- spline(predPhi$Time, predPhi$ucl, n=100)$y
p <- ggplot(aes(x=x, y=y), data=pred.sm)
p <- p + geom_line() + geom_ribbon(aes(ymin=lcl, ymax=ucl), alpha=0.25) +
  ggtitle("Female dipper survival\n") + xlab("Time") + ylab("Phi(t)")
p


The default is 3 degrees of freedom and the above will work with the existing version of RMark. However, if you want to specify the degrees of freedom as in ~bs(Time,df=4) then you need version 2.1.5 of RMark at:
https://docs.google.com/folder/d/0B77g1ScdUwVeOVJNUVVGS0YtWE0/edit?usp=sharing
I had to make the output extraction code smarter so it didn't trip on the = sign. You need the splines package to use bs and for the above code you need ggplot2 to produce the plot.

I keep the newest version of RMark and other software at the above location. To use the new version, you need to download the appropriate zip file and then use in the R menu, Packages/Install package from local zip. Browse to the location of the downloaded file and select it to install.

regards --jeff

[jlaake]

I meant to add the following to my previous post. This approach can be used with the MARK interface by using the splines package in R with RMark to create the design matrix and then it can be copied into the DM in the MARK interface as long as the PIM structure you specify in MARK matches what you use with RMark. Obviously this is easier in RMark but now Evan can scratch this off as a perceived advantage of RMark over MARK interface.

Code: Select all

> library(splines)
> library(RMark)
> data(dipper)
> dipper.proc=process.data(dipper,model="CJS")
> dipper.ddl=make.design.data(dipper.proc)
> model.matrix(~bs(Time),dipper.ddl$Phi)
   (Intercept) bs(Time)1 bs(Time)2 bs(Time)3
1            1     0.000     0.000     0.000
2            1     0.384     0.096     0.008
3            1     0.432     0.288     0.064
4            1     0.288     0.432     0.216
5            1     0.096     0.384     0.512
6            1     0.000     0.000     1.000
7            1     0.384     0.096     0.008
8            1     0.432     0.288     0.064
9            1     0.288     0.432     0.216
10           1     0.096     0.384     0.512
11           1     0.000     0.000     1.000
12           1     0.432     0.288     0.064
13           1     0.288     0.432     0.216
14           1     0.096     0.384     0.512
15           1     0.000     0.000     1.000
16           1     0.288     0.432     0.216
17           1     0.096     0.384     0.512
18           1     0.000     0.000     1.000
19           1     0.096     0.384     0.512
20           1     0.000     0.000     1.000
21           1     0.000     0.000     1.000
attr(,"assign")
[1] 0 1 1 1
>