predict_real and mlogit parameters

[jlaake]

In working offlist with a person, I discovered that predict_real does not work with mlogit parameters. Not sure when I'll get around to fixing it (if ever.- getting old). But it bugged me that it didn't work so I just devised a way to get predictions from covariate.predictions with specified design data covariate values.  This will work with non mlogit parameters as well. Had I thought of this before I wouldn't have written predict_real. You do have to know some R to modify but it does work. You do have to be careful if the covariate (cov here) is used in an interaction because you are modifying values in the design matrix. But for straightforward use of a design covariate this will work. It should be possible to feed these estimates to TransitionMatrix if you want estimates for the subtracted stratum (typically probability of staying for Psi).

Code: Select all

# use mstrata data and fit a model with design covariate cov in Psi
# to create an example to work with.
library(RMark)
data(mstrata)
dp=process.data(mstrata,model="Multistrata")
ddl=make.design.data(dp)
ddl$Psi$cov=1:36
model=mark(dp,ddl,model.parameters=list(Psi=list(formula=~stratum+cov)))
# save model into model.sub
model.sub=model
# This is a loop for a set of values to plug into cov in the design matrix and get real predictions for Psi
# and save in a list (est.list) 
# specify set of values 
values=1:36
# create empty list - one for each of the values
est.list=list(length=length(values))
# loop over values
for(i in 1:length(values))
{ 
   # get original design matrix from the model
   dm=model$design.matrix 
   # find the column for cov in Psi
   icol=which(colnames(dm)=="Psi:cov") 
   # find non-zero rows
   irows=which(dm[,icol]!="0") 
   # substitute the value into the design matrix
   dm[irows,icol]=values[i] 
   # put the new design matrix into the model.sub design matrix
   model.sub$design.matrix=dm
   # call covariate.predictions and store the estimates with this value
   est.list[[i]]=covariate.predictions(model.sub,indices=37:72)$estimates
   i=i+1
}