Multi-state transition SEs (delta method?)

[awan]

Dear all,

We have a multi-state live dead model that estimates survival and transitions among three different breeding states. We also assess how these are affected by environmental covariates and use covariate.predictions to plot these relationships. When calling the indices for the covariate.predictions, it is obviously only possible to include transitions among others states, whilst the probability of remaining in the same state is 1 - x1 - x2. We have been asked to included standard errors for remaining in the same state, and a suggestion was to use the delta method. We have been reading Appendix B, the workshop notes and various topics on the forum but the context tends to be different (e.g. most focus is on product of estimates). I thought that perhaps our example is similar to applying the delta method with "sum", where it is actually a subtraction. So for example, following this example, https://www.montana.edu/rotella/documen ... 4RMark.rmd, I tried

Code: Select all

rr <- get.real(ms.results[[6]], "Psi", se = TRUE, vcv = TRUE) ##reals of top model
sigma <- rr$vcv.real
SN2002 <- rr$estimates$estimate[1] ##first time stamp for SN
SF2002 <- rr$estimates$estimate[121] ##first tim stamp for SF
Diff <- 1 - SN2002 - SF2002
seDiff <- deltamethod(~1 - x1 - x2, c(SN2002, SF2002),sigma)


Here I meet a couple of issues. One is that we have 15 years of data, and three states, meaning it is a quite a large matrix. Meanwhile, the delta method function is expecting a 2x2 matrix. So I retried the above again, but this time created my own matrix using the appropriate values from the vcv as follows

Code: Select all

rr <- get.real(ms.results[[6]], "Psi", se = TRUE, vcv = TRUE)
sigma <- matrix(c(rr$vcv.real[1,1], rr$vcv.real[121,1], rr$vcv.real[121,1], rr$vcv.real[1,1]), nrow = 2, ncol = 2) ##vcv matrix is thus r1 = 0.0005912684, -0.0003690198, r2 = -0.0003690198, 0.0005912684
##note that rr$vcv.real[121,1] = rr$vcv.real[1,121]
SN2002 <- rr$estimates$estimate[1]  ##S to N for 2002, value is 0.2339315
SF2002 <- rr$estimates$estimate[121] ##S to F for 2002, value is 0.1520598
Diff <- 1 - SN2002 - SF2002 ##Mlogit hence S to S is 1 - SF - SN
seDiff <- deltamethod(~1 - x1 - x2, c(SN2002, SF2002),sigma)


The above gives an SE of 0.0210831, which is the region of the SE of SN (0.02431601) and SF (0.02725081). However, if I change the above formula to simply x1 - x2, the SE is 0.04382438 - which happens to be the same as the comparison on the above link, i.e.

Code: Select all

(se.diff.direct=sqrt(sum(diag(sigma))-2*sigma[1,2])).


Lastly, the SE of deltamethod(~1 - x1 - x2) = the SE of deltamethod(~x1 + x2).

Firstly, is my current notation correct for implementing and specifying the formula of the delta method. Secondly, will I need to rinse and repeat this process 45 times (i.e. 3 states x 15 years), creating a new vcv matrix everytime to obtain these SEs? Is there a more efficient approach that I am missing?

Many thanks in advance for any guidance you are able to provide.

[jlaake]

Look at help for function TransitionMatrix. It should do what you want.

[awan]

Thank you Jeff for the prompt response. That is indeed a much quicker and easier way of getting the SEs for each transitions (including staying in the same state).

Of course, upon completion of this step, I realise I hit a stumbling block in that this doesn't help me get towards the SE of what would effectively be a beta estimate (in order to effectively plot the equivalent of the covariate predictions for remaining in the same state, i.e. the relationship between the transition and environmental covariate). Do you have any ideas or if you can point me in the right direction I would be grateful.

p.s. we have three environmental covariates, so when using covariate.predictions to plot the relationships, we do it one covariate at a time by setting the other two covariates to zero (mean is zero). This creates a challenge for using the reals because it means the transitions are based on three covariates. To use covariate.predictions we have obviously added the covariates as time-varying individual covariates.

[jlaake]

Sorry for the delay. One of the issues (benefits) of retirement is that there are a lot more demands on my time other than work. Anyhow, I have a partial solution for you below using the mstrata example in RMark. I say partial because I only use a single covariate. You will have to modify the code to handle 3 covariates. I suggest that you look at "fc" for your model and you'll see it has multiple entries for each of the covariates, so you'll have to specify which rows you want to set for each covariate and the value you want to assign. You'll also have to modify the code for "df" to include a column for each covariate. But other than that this should get you where you want to go. But before you start to modify it, try with this example so you can understand what it is doing. mstrata doesn't have covariates, so I just created a dummy one.

Code: Select all

data(mstrata)
mstrata$cov=rnorm(nrow(mstrata),0,3)
run.mstrata=function()
{
#
# Process data
#
mstrata.processed=process.data(mstrata,model="Multistrata")
#
# Create default design data
#
mstrata.ddl=make.design.data(mstrata.processed)
#
#  Define range of models for S; note that the betas will differ from the output
#  in MARK for the ~stratum = S(s) because the design matrix is defined using
#  treatment contrasts for factors so the intercept is stratum A and the other
#  two estimates represent the amount that survival for B abd C differ from A.
#  You can use force the approach used in MARK with the formula ~-1+stratum which
#  creates 3 separate Betas - one for A,B and C.
#
S.stratum=list(formula=~stratum)
#
#  Define range of models for p
#
p.stratum=list(formula=~stratum)
#
#  Define range of models for Psi; what is denoted as s for Psi
#  in the Mark example for Psi is accomplished by -1+stratum:tostratum which
#  nests tostratum within stratum.  Likewise, to get s*t as noted in MARK you
#  want ~-1+stratum:tostratum:time with time nested in tostratum nested in
#  stratum.
#
Psi.s=list(formula=~-1+stratum:tostratum+cov)
#
# Create model list and run assortment of models
#
model.list=create.model.list("Multistrata")
#
# Run the list of models
#
mstrata.results=mark.wrapper(model.list,data=mstrata.processed,ddl=mstrata.ddl)
#
# Return model table and list of models
#
return(mstrata.results)
}
mstrata.results=run.mstrata()
mstrata.results
model=mstrata.results[[1]]

# find entries in design matrix that are individual covariates
fc=find.covariates(model)

# loop over different values of covariate and compute reals and transition values
# and store in dataframe
df=NULL
for(v in -3:3)
{
  fc$value=v
  dm=fill.covariates(model,fc)
  Psilist=get.real(model,"Psi",vcv=TRUE,design=dm)
  Psivalues=Psilist$estimates
  TransitionMatrix(Psivalues[Psivalues$time==1,])
# call it again but specify the vc matrix to get se and conf interval
  tm=TransitionMatrix(Psivalues[Psivalues$time==1,],vcv.real=Psilist$vcv.real)
  from=rep(rownames(tm$TransitionMat),ncol(tm$TransitionMat))
  to=rep(rownames(tm$TransitionMat),each=ncol(tm$TransitionMat))
  df=rbind(df,data.frame(cov=v,from=from,to=to,transProb=as.vector(tm$TransitionMat),transProb.se=as.vector(tm$se.TransitionMat),
                transProb.lcl=as.vector(tm$lcl.TransitionMat),transProb.ucl=as.vector(tm$ucl.TransitionMat)))
}
df



[awan]

Hi Jeff,

Of course no problem, and I appreciate how you continue to provide your valued support even during your retirement!!

In the meantime, I had also been thinking a bit more about the problem, and realised that the results of covariate.predictions also includes a vcv matrix and hence it may be possible to apply the delta method to the predictions (where a range of values is supplied for one covariate whilst the other two are held at the mean, in this case zero).

I am sure my code is not very efficient, especially as it is modified slightly for each covariate, and for the three different transitions of SS, FF and NN (and hence needs to call on different line numbers in the results object), but I thought I would share what I have done. I will also try have a look at your shared example to see how the estimates compare. In the meantime, this is the code I used to implement the delta method on the predictions:

Code: Select all

##indices are six transitions.
##Code could probably be streamlined by repeating per S->, N-> and F-> indices
##Here I do temperature only, need to be repeated for other two covariates too.
TempPred <- covariate.predictions(ms.TrantoPlot[[1]], data = data.frame(temp2002 = seq(-2, 2, 0.1)),
indices = c(769,784,799,814,829,844), drop = FALSE)

##Extract real estimates and VCV matrix
ToPlotTemp <- TempPred$estimates
ToPlotTempVCV <- TempPred$vcv
ToPlotTemp$Covariate <- "Temperature"
ToPlotTemp <- ToPlotTemp[order(ToPlotTemp$par.index),]

##create object to store delta results
TempDelta <- data.frame(Transition = rep(c("StoS", "NtoN", "FtoF"),each = 41), NewSE = NA)
library(msm)
##estimate transitions for SS (note that rows 1:41 is S->N and 42:82 is S->F)
for(i in 1:41){
##to create matrix,
##use the info in the vcv.index column to call the correct row/column in VCV object
sigma <- matrix(c( 
ToPlotTempVCV[ToPlotTemp$vcv.index[i], ToPlotTemp$vcv.index[i]],
ToPlotTempVCV[ToPlotTemp$vcv.index[i], ToPlotTemp$vcv.index[41+i]],
ToPlotTempVCV[ToPlotTemp$vcv.index[41+i], ToPlotTemp$vcv.index[i]],
ToPlotTempVCV[ToPlotTemp$vcv.index[41+i], ToPlotTemp$vcv.index[41+i]]), nrow = 2, ncol = 2)
SN <- ToPlotTemp$estimate[i]  ##S to N estimate
SF <- ToPlotTemp$estimate[41+i] ##S to F estimate
Diff <- 1 - SN - SF ##Mlogit hence S to S is 1 - SF - SN
seDiff <- deltamethod(~1 - x1 - x2, c(SN2002, SF2002),sigma)
TempDelta$NewSE[i] <- seDiff
}
##Then create new rows in TempPred object for SS estimates
##Paste in the TempDelta$NewSE
##Note that lcl and ucl are now uniform (i.e. estimate +- 1.96 * se)
##Normally that would not be the case when converting from the betas


This is the end results of the image I was trying to create, where S to S now includes CIs.

[jlaake]

You can do it that way as well but I think if you look at my code it would end up being more efficient as it computes all of the transitions and their std errors and confidence intervals.

But let me ask a different question. Why are you using time-varying individual covariates? Is temperature varying by individual? If not then temperature should be used as a design covariate. Something I added in 2016 v2.2.0 was predict_real which allows you to specify a range of values for design covariates. Now it is not necessary for you to change because what you are doing should work. But individual covariates increase model fitting times in comparison to design covariates.

--jeff

[jlaake]

If it helps, here is an example with 3 covariates. In each case, I vary one covariate at a time with the other two held at 0 value. But you can nest the loops and compute all of the sets of values - in this case a 3-d cube over -3 to 3 for each variable. You would obviously use different ranges relevant to your variables and different variable names. Displaying the results will be challenging but you could use 3-d plots with pairs of covariates.

Code: Select all

data(mstrata)
mstrata$cov1=rnorm(nrow(mstrata),0,3)
mstrata$cov2=rnorm(nrow(mstrata),0,3)
mstrata$cov3=rnorm(nrow(mstrata),0,3)
run.mstrata=function()
{
#
# Process data
#
mstrata.processed=process.data(mstrata,model="Multistrata")
#
# Create default design data
#
mstrata.ddl=make.design.data(mstrata.processed)
#
#  Define range of models for S; note that the betas will differ from the output
#  in MARK for the ~stratum = S(s) because the design matrix is defined using
#  treatment contrasts for factors so the intercept is stratum A and the other
#  two estimates represent the amount that survival for B abd C differ from A.
#  You can use force the approach used in MARK with the formula ~-1+stratum which
#  creates 3 separate Betas - one for A,B and C.
#
S.stratum=list(formula=~stratum)
#
#  Define range of models for p
#
p.stratum=list(formula=~stratum)
#
#  Define range of models for Psi; what is denoted as s for Psi
#  in the Mark example for Psi is accomplished by -1+stratum:tostratum which
#  nests tostratum within stratum.  Likewise, to get s*t as noted in MARK you
#  want ~-1+stratum:tostratum:time with time nested in tostratum nested in
#  stratum.
#
Psi.s=list(formula=~-1+stratum:tostratum+cov1+cov2+cov3)
#
# Create model list and run assortment of models
#
model.list=create.model.list("Multistrata")
#
# Run the list of models
#
mstrata.results=mark.wrapper(model.list,data=mstrata.processed,ddl=mstrata.ddl)
#
# Return model table and list of models
#
return(mstrata.results)
}
mstrata.results=run.mstrata()
mstrata.results
model=mstrata.results[[1]]

# find entries in design matrix that are individual covariates
fc=find.covariates(model)

# loop over different values of covariate cov1 holding cov2=0 and compute reals and transition values
# and store in dataframe
df=NULL
for(v in -3:3)
{
  fc$value[fc$var=="cov1"]=v
  fc$value[fc$var=="cov2"]=0
  fc$value[fc$var=="cov3"]=0
  dm=fill.covariates(model,fc)
  Psilist=get.real(model,"Psi",vcv=TRUE,design=dm)
  Psivalues=Psilist$estimates
  TransitionMatrix(Psivalues[Psivalues$time==1,])
# call it again but specify the vc matrix to get se and conf interval
  tm=TransitionMatrix(Psivalues[Psivalues$time==1,],vcv.real=Psilist$vcv.real)
  from=rep(rownames(tm$TransitionMat),ncol(tm$TransitionMat))
  to=rep(rownames(tm$TransitionMat),each=ncol(tm$TransitionMat))
  df=rbind(df,data.frame(cov1=v,cov2=0,cov3=0,from=from,to=to,transProb=as.vector(tm$TransitionMat),transProb.se=as.vector(tm$se.TransitionMat),
                transProb.lcl=as.vector(tm$lcl.TransitionMat),transProb.ucl=as.vector(tm$ucl.TransitionMat)))
}

# loop over different values of covariate cov1 holding cov2=0 and compute reals and transition values
# and store in dataframe
for(v in -3:3)
{
  fc$value[fc$var=="cov1"]=0
  fc$value[fc$var=="cov2"]=v
  fc$value[fc$var=="cov3"]=0
  dm=fill.covariates(model,fc)
  Psilist=get.real(model,"Psi",vcv=TRUE,design=dm)
  Psivalues=Psilist$estimates
  TransitionMatrix(Psivalues[Psivalues$time==1,])
  # call it again but specify the vc matrix to get se and conf interval
  tm=TransitionMatrix(Psivalues[Psivalues$time==1,],vcv.real=Psilist$vcv.real)
  from=rep(rownames(tm$TransitionMat),ncol(tm$TransitionMat))
  to=rep(rownames(tm$TransitionMat),each=ncol(tm$TransitionMat))
  df=rbind(df,data.frame(cov1=0,cov2=v,cov3=0,from=from,to=to,transProb=as.vector(tm$TransitionMat),transProb.se=as.vector(tm$se.TransitionMat),
                         transProb.lcl=as.vector(tm$lcl.TransitionMat),transProb.ucl=as.vector(tm$ucl.TransitionMat)))
}


for(v in -3:3)
{
  fc$value[fc$var=="cov1"]=0
  fc$value[fc$var=="cov2"]=0
  fc$value[fc$var=="cov3"]=v
  dm=fill.covariates(model,fc)
  Psilist=get.real(model,"Psi",vcv=TRUE,design=dm)
  Psivalues=Psilist$estimates
  TransitionMatrix(Psivalues[Psivalues$time==1,])
  # call it again but specify the vc matrix to get se and conf interval
  tm=TransitionMatrix(Psivalues[Psivalues$time==1,],vcv.real=Psilist$vcv.real)
  from=rep(rownames(tm$TransitionMat),ncol(tm$TransitionMat))
  to=rep(rownames(tm$TransitionMat),each=ncol(tm$TransitionMat))
  df=rbind(df,data.frame(cov1=0,cov2=0,cov3=v,from=from,to=to,transProb=as.vector(tm$TransitionMat),transProb.se=as.vector(tm$se.TransitionMat),
                         transProb.lcl=as.vector(tm$lcl.TransitionMat),transProb.ucl=as.vector(tm$ucl.TransitionMat)))
}

df


[awan]

Hi Jeff,

Thank you for the updated example using three covariates. I will try applying it to our data.

In response to your previous post, I had missed the update whereby one could start using predict_real instead. I had indeed originally included it as a design covariate, but my "current knowledge" was based on the previous process of including the covariates as time-varying individual covariates so that the function covariate.predictions could be used. This update will make my mark data file much simpler I guess this is one of the challenges of google searches or searching the forum. Doing a search on plotting environmental covariates and I hit this topic viewtopic.php?f=21&t=3044&p=9773&hilit=time+varying+covariates#p9773 and hence missed the subsequent updates to RMark.