prep.regimes working?

[mhunter1]

Original report by Sy-Miin Chow (Bitbucket: symiin, GitHub: symiin).

---

Prep.regimes doesn't seem to be working correctly with covariates.

I specified this model where the covariates x1 and x2 are used to predict the LO parameters associated with the (1,1) and (2,1) elements of the transition probability matrix. When I printed out the specified model, I got
Slot "covariates":
character(0)

When optimizing the full model (RS-linearODE), the parameters for the covariates don't seem to be updated at all. Looking at prep.regimes, I also didn't see how the covariates were actually used inside prep.regimes.

#!R
# Regime-switching function
# The RS model assumes that each element of the transition probability 
# matrix (TPM) can be expressed as a linear predictor (lp).
# LPM = 
# lp(p11) ~ 1 + x1 + x2 + ... + xn,   lp(p12) ~ 1 + x1 + x2 + ... + xn
# lp(p21) ~ 1 + x1 + x2 + ... + xn,   lp(p22) ~ 1 + x1 + x2 + ... + xn
# Here I am specifying lp(p11) and lp(p21); the remaining elements
# lp(p12) and lp(p22) are fixed at zero.

regimes <- prep.regimes(
  values=matrix(c(6,.5,-.3,rep(0,3),
                  -3,-1.5,-1,rep(0,3)), 
                nrow=2, ncol=6,byrow=T), # nrow=numRegimes, ncol=numRegimes*(numCovariates+1)
  params=matrix(c("a_{11}","d_{11,1}","d_{11,2}",rep("fixed",3),
                  "a_{21}","d_{21,1}","d_{21,2}",rep("fixed",3),), 
                nrow=2, ncol=6,byrow=T), covariates=c('x1', 'x2'))

[mhunter1]

Original comment by Michael Hunter (Bitbucket: mhunter, GitHub: mhunter).

---

There is at least one typo in your code. There's an extra comma in the params matrix. Here's a corrected version.

#!R
regimes <- prep.regimes(
  values=matrix(c(6,.5,-.3,rep(0,3),
                  -3,-1.5,-1,rep(0,3)), 
                nrow=2, ncol=6,byrow=T), # nrow=numRegimes, ncol=numRegimes*(numCovariates+1)
  params=matrix(c("a_{11}","d_{11,1}","d_{11,2}",rep("fixed",3),
                  "a_{21}","d_{21,1}","d_{21,2}",rep("fixed",3)), 
                nrow=2, ncol=6,byrow=T), covariates=c('x1', 'x2'))

The covariates slot is not currently populated. Rather, it is assumed that the C co_variates argument common to many of these functions is the same one referred to in prep.regimes. I believe issue #34 is intended to update this.

[mhunter1]

Original comment by Sy-Miin Chow (Bitbucket: symiin, GitHub: symiin).

---

Fixing the typo changed the locations of the parameters that didn't get updated in the optimization, but still not working (some parameters still didn't get updated). Do you see any other errors in my use of prep.regimes? My covariate names already matched the names of the covariates I sent to dynr.data.

colnames(thedata) = c("ID","Time","y1","y2","x1","x2")
data <- dynr.data(thedata, id="ID", time="Time",observed=paste0('y', 1:2),
covariates=paste0('x', 1:2))

Also, doing cat(writeCcode(regimes)$c.string) gave error:
Error in params[i] - 1 : non-numeric argument to binary operator

[mhunter1]

Original comment by Michael Hunter (Bitbucket: mhunter, GitHub: mhunter).

---

I don't see anything obvious. Maybe we're having trouble with your parameter names? The underscores, commas, and the curly braces might be throwing something off.

You should be able to inspect the C code for the regime-switching with

#!R
cat(model$regimes$c.string)

provided that model was created by dynr.model

Try

#!R
regimes <- prep.regimes(
  values=matrix(c(6,.5,-.3,rep(0,3),
                  -3,-1.5,-1,rep(0,3)), 
                nrow=2, ncol=6,byrow=T), # nrow=numRegimes, ncol=numRegimes*(numCovariates+1)
  params=matrix(c("a11","d111","d112",rep("fixed",3),
                  "a21","d211","d212",rep("fixed",3)), 
                nrow=2, ncol=6,byrow=T), covariates=c('x1', 'x2'))

[mhunter1]

Original comment by Sy-Miin Chow (Bitbucket: symiin, GitHub: symiin).

---

The generated C code doesn't use the covariates correctly. See below. Can you check your function for converting from prep.regimes to c code to make sure that the covariates are being used correctly?

cat(model$regimes$c.string)
void function_regime_switch(size_t t, size_t type, double *param, const gsl_vector *co_variate, gsl_matrix *regime_switch_mat){
gsl_vector *Pvector = gsl_vector_calloc(2);

gsl_vector *Presult = gsl_vector_calloc(2);

gsl_vector_set(Pvector, 0, param[7]);
gsl_vector_set(Pvector, 1, param[9]);

mathfunction_softmax(Pvector, Presult);
gsl_matrix_set_row(regime_switch_mat, 0, Presult);

gsl_vector_set_zero(Pvector);
gsl_vector_set(Pvector, 0, param[11]);

mathfunction_softmax(Pvector, Presult);
gsl_matrix_set_row(regime_switch_mat, 1, Presult);

gsl_vector_set_zero(Pvector);
gsl_vector_free(Pvector);

gsl_vector_free(Presult);

}

[mhunter1]

Original comment by Michael Hunter (Bitbucket: mhunter, GitHub: mhunter).

---

Argh! There was a very tiny change that now produces

#!R
cat(model$regimes$c.string)

#!c
void function_regime_switch(size_t t, size_t type, double *param, const gsl_vector *co_variate, gsl_matrix *regime_switch_mat){
        gsl_matrix *Gmatrix = gsl_matrix_calloc(2, 2);

        gsl_vector *Pvector = gsl_vector_calloc(2);

        gsl_vector *Presult = gsl_vector_calloc(2);

        gsl_vector_set(Pvector, 0, param[7]);

        gsl_matrix_set(Gmatrix, 0, 0, param[8]);
        gsl_matrix_set(Gmatrix, 0, 1, param[9]);

        gsl_blas_dgemv(CblasNoTrans, 1.0, Gmatrix, co_variate, 1.0, Pvector);

        mathfunction_softmax(Pvector, Presult);
        gsl_matrix_set_row(regime_switch_mat, 0, Presult);

        gsl_vector_set_zero(Pvector);
        gsl_matrix_set_zero(Gmatrix);
        gsl_vector_set(Pvector, 1, param[10]);

        gsl_matrix_set(Gmatrix, 1, 0, param[11]);
        gsl_matrix_set(Gmatrix, 1, 1, param[12]);

        gsl_blas_dgemv(CblasNoTrans, 1.0, Gmatrix, co_variate, 1.0, Pvector);

        mathfunction_softmax(Pvector, Presult);
        gsl_matrix_set_row(regime_switch_mat, 1, Presult);

        gsl_vector_set_zero(Pvector);
        gsl_matrix_set_zero(Gmatrix);
        gsl_matrix_free(Gmatrix);

        gsl_vector_free(Pvector);

        gsl_vector_free(Presult);

}

Does there version updated with 84a963f work for you? The version in the repo still has to have the transformation function updated.

[mhunter1]

Original comment by Sy-Miin Chow (Bitbucket: symiin, GitHub: symiin).

---

Yup, with the change now RSlinearODe.r works. I added some error checking to prep.regimes and updated the roxygen slightly. But need Sukruth to beautify.