TABLE OF CONTENTS

CMetropolisHastings/MH_Regression [ Functions ]

NAME

    MH_Regression --- MH step for regression coefficients

FUNCTION

Update the regression coefficients by Metropolis-Hastings. See also mh.coef.

SYNOPSIS

2012 void MH_Regression(curveP hazard, curveP frailty, regressionP regression)

INPUTS

    hazard        CCurve for the hazard
    frailty       CCurve for the frailty
    regression    CRegression structure

SOURCE

2016 {
2017     double baselik, candlik;
2018     // base likelihood
2019     baselik = LikelihoodRegression(hazard, frailty, regression);
2020     double * cand = (double *) calloc( regression->p, sizeof(double));
2021     double * oldlp = (double *) malloc( regression->n * sizeof(double));
2022     double * oldelp = (double *) malloc( regression->n * sizeof(double));
2023     double * oldfrailelp = (double *) malloc( regression->n * sizeof(double));
2024     double * oldcoef = (double *) malloc( regression->p * sizeof(double));
2025     // store the old regression information
2026     char trans = 'N'; double c0 = 0; int c1 = 1; double c1d = 1;
2027     dcopyWrapper(regression->p, regression->coefficients, oldcoef);
2028     dcopyWrapper(regression->n, regression->lp, oldlp);
2029     dcopyWrapper(regression->n, regression->elp, oldelp);
2030     dcopyWrapper(regression->n, regression->frailelp, oldfrailelp);
2031     //generate candidate parameters
2032     mvrnorm(regression->p, cand, regression->coefficients, regression->CholCov, regression->tun[0]); 
2033     // Change the regression object with the new lp and elps
2034     dcopyWrapper(regression->p, cand, regression->coefficients);
2035     F77_CALL(dgemv)(&trans, &(regression->n), &(regression->p), &c1d, regression->covariates,
2036             &(regression->n), regression->coefficients, &c1, &c0, regression->lp, &c1);
2037     for(int i=0; i < regression->n; i++) regression->elp[i] = exp(regression->lp[i]);
2038     diagmvWrapper(regression->n, regression->frailrep, regression->elp, regression->frailelp);
2039     candlik = LikelihoodRegression(hazard, frailty, regression);
2040     int acc = AcceptReject(baselik, candlik, 1);
2041     if(!acc){
2042         // not accepted, so restore old regression
2043         dcopyWrapper(regression->p, oldcoef, regression->coefficients);
2044         dcopyWrapper(regression->n, oldlp, regression->lp);
2045         dcopyWrapper(regression->n, oldelp, regression->elp);
2046         dcopyWrapper(regression->n, oldfrailelp, regression->frailelp);
2047     }
2048     regression->Accept[0] = (double) acc; 
2049     free(cand);
2050     free(oldlp);
2051     free(oldelp);
2052     free(oldfrailelp);
2053     free(oldcoef);
2054 }