version 1.225, 2016/07/12 08:40:03
|
version 1.226, 2016/07/12 18:42:34
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
$State$ |
$State$ |
$Log$ |
$Log$ |
|
Revision 1.226 2016/07/12 18:42:34 brouard |
|
Summary: temp |
|
|
Revision 1.225 2016/07/12 08:40:03 brouard |
Revision 1.225 2016/07/12 08:40:03 brouard |
Summary: saving but not running |
Summary: saving but not running |
|
|
Line 689 Back prevalence and projections:
|
Line 692 Back prevalence and projections:
|
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying |
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying |
nhstepm*hstepm matrices. Returns p3mat[i][j][h] after calling |
nhstepm*hstepm matrices. Returns p3mat[i][j][h] after calling |
p3mat[i][j][h]=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ |
p3mat[i][j][h]=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ |
1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); |
1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); |
|
|
|
Important routines |
|
|
|
- func (or funcone), computes logit (pij) distinguishing |
|
o fixed variables (single or product dummies or quantitative); |
|
o varying variables by: |
|
(1) wave (single, product dummies, quantitative), |
|
(2) by age (can be month) age (done), age*age (done), age*Vn where Vn can be: |
|
% fixed dummy (treated) or quantitative (not done because time-consuming); |
|
% varying dummy (not done) or quantitative (not done); |
|
- Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities) |
|
and returns the number of efficient covariates cptcoveff and modalities nbcode[Tvar[k]][1]= 0 and nbcode[Tvar[k]][2]= 1 usually. |
|
- printinghtml which outputs results like life expectancy in and from a state for a combination of modalities of dummy variables |
|
o There are 2*cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, éliminating 1 1 if |
|
race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless. |
|
|
|
|
|
|
Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). |
Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). |
Institut national d'études démographiques, Paris. |
Institut national d'études démographiques, Paris. |
This software have been partly granted by Euro-REVES, a concerted action |
This software have been partly granted by Euro-REVES, a concerted action |
Line 761 Back prevalence and projections:
|
Line 781 Back prevalence and projections:
|
#include <stdio.h> |
#include <stdio.h> |
#include <stdlib.h> |
#include <stdlib.h> |
#include <string.h> |
#include <string.h> |
|
#include <ctype.h> |
|
|
#ifdef _WIN32 |
#ifdef _WIN32 |
#include <io.h> |
#include <io.h> |
Line 1015 double ***cotvar; /* Time varying covari
|
Line 1036 double ***cotvar; /* Time varying covari
|
double ***cotqvar; /* Time varying quantitative covariate itqv */ |
double ***cotqvar; /* Time varying quantitative covariate itqv */ |
double idx; |
double idx; |
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ |
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ |
int *Typevar; /**< 1 for qualitative fixed, 2 for quantitative fixed, 3 for qualitive varying, 4 for quanti varying*/ |
int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */ |
|
int *Fixed; /** Fixed[Tvar[k]] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */ |
|
int *Dummy; /** Dummy[Tvar[k]] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */ |
int *Tage; |
int *Tage; |
int *Ndum; /** Freq of modality (tricode */ |
int *Ndum; /** Freq of modality (tricode */ |
/* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */ |
/* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */ |
Line 2890 double ***hbxij(double ***po, int nhstep
|
Line 2913 double ***hbxij(double ***po, int nhstep
|
/*************** log-likelihood *************/ |
/*************** log-likelihood *************/ |
double func( double *x) |
double func( double *x) |
{ |
{ |
int i, ii, j, k, mi, d, kk; |
int i, ii, j, k, mi, d, kk; |
int ioffset=0; |
int ioffset=0; |
double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; |
double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1]; |
double **out; |
double **out; |
double sw; /* Sum of weights */ |
double sw; /* Sum of weights */ |
double lli; /* Individual log likelihood */ |
double lli; /* Individual log likelihood */ |
int s1, s2; |
int s1, s2; |
int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quatitative time varying covariate */ |
int iv=0, iqv=0, itv=0, iqtv=0 ; /* Index of varying covariate, fixed quantitative cov, time varying covariate, quatitative time varying covariate */ |
double bbh, survp; |
double bbh, survp; |
long ipmx; |
long ipmx; |
double agexact; |
double agexact; |
/*extern weight */ |
/*extern weight */ |
/* We are differentiating ll according to initial status */ |
/* We are differentiating ll according to initial status */ |
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ |
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ |
/*for(i=1;i<imx;i++) |
/*for(i=1;i<imx;i++) |
printf(" %d\n",s[4][i]); |
printf(" %d\n",s[4][i]); |
*/ |
*/ |
|
|
++countcallfunc; |
++countcallfunc; |
|
|
cov[1]=1.; |
cov[1]=1.; |
|
|
for(k=1; k<=nlstate; k++) ll[k]=0.; |
for(k=1; k<=nlstate; k++) ll[k]=0.; |
ioffset=0; |
ioffset=0; |
if(mle==1){ |
if(mle==1){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
/* Computes the values of the ncovmodel covariates of the model |
/* Computes the values of the ncovmodel covariates of the model |
depending if the covariates are fixed or varying (age dependent) and stores them in cov[] |
depending if the covariates are fixed or varying (age dependent) and stores them in cov[] |
Then computes with function pmij which return a matrix p[i][j] giving the elementary probability |
Then computes with function pmij which return a matrix p[i][j] giving the elementary probability |
to be observed in j being in i according to the model. |
to be observed in j being in i according to the model. |
*/ |
*/ |
ioffset=2+nagesqr+cptcovage; |
ioffset=2+nagesqr+cptcovage; |
/* for (k=1; k<=cptcovn;k++){ /\* Simple and product covariates without age* products *\/ */ |
/* for (k=1; k<=cptcovn;k++){ /\* Simple and product covariates without age* products *\/ */ |
for (k=1; k<=ncoveff;k++){ /* Simple and product covariates without age* products */ |
for (k=1; k<=ncoveff;k++){ /* Simple and product covariates without age* products */ |
cov[++ioffset]=covar[Tvar[k]][i]; |
cov[++ioffset]=covar[Tvar[k]][i]; |
} |
} |
for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitatives and Fixed covariates */ |
for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitatives and Fixed covariates */ |
cov[++ioffset]=coqvar[iqv][i]; |
cov[++ioffset]=coqvar[iqv][i]; |
} |
} |
|
|
/* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] |
/* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] |
is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] |
is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] |
has been calculated etc */ |
has been calculated etc */ |
/* For an individual i, wav[i] gives the number of effective waves */ |
/* For an individual i, wav[i] gives the number of effective waves */ |
/* We compute the contribution to Likelihood of each effective transition |
/* We compute the contribution to Likelihood of each effective transition |
mw[mi][i] is real wave of the mi th effectve wave */ |
mw[mi][i] is real wave of the mi th effectve wave */ |
/* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i]; |
/* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i]; |
s2=s[mw[mi+1][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i] |
And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i] |
But if the variable is not in the model TTvar[iv] is the real variable effective in the model: |
But if the variable is not in the model TTvar[iv] is the real variable effective in the model: |
meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i] |
meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i] |
*/ |
*/ |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */ |
for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */ |
cov[ioffset+itv]=cotvar[mw[mi][i]][itv][i]; |
cov[ioffset+itv]=cotvar[mw[mi][i]][itv][i]; |
} |
} |
for(iqtv=1; iqtv <= nqtveff; iqtv++){ /* Varying quantitatives covariates */ |
for(iqtv=1; iqtv <= nqtveff; iqtv++){ /* Varying quantitatives covariates */ |
if(cotqvar[mw[mi][i]][iqtv][i] == -1){ |
if(cotqvar[mw[mi][i]][iqtv][i] == -1){ |
printf("i=%d, mi=%d, iqtv=%d, cotqvar[mw[mi][i]][iqtv][i]=%f",i,mi,iqtv,cotqvar[mw[mi][i]][iqtv][i]); |
printf("i=%d, mi=%d, iqtv=%d, cotqvar[mw[mi][i]][iqtv][i]=%f",i,mi,iqtv,cotqvar[mw[mi][i]][iqtv][i]); |
} |
} |
cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][iqtv][i]; |
cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][iqtv][i]; |
} |
} |
/* ioffset=2+nagesqr+cptcovn+nqv+ntv+nqtv; */ |
/* ioffset=2+nagesqr+cptcovn+nqv+ntv+nqtv; */ |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
} |
} |
for(d=0; d<dh[mi][i]; d++){ |
for(d=0; d<dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
cov[2]=agexact; |
cov[2]=agexact; |
if(nagesqr==1) |
if(nagesqr==1) |
cov[3]= agexact*agexact; /* Should be changed here */ |
cov[3]= agexact*agexact; /* Should be changed here */ |
for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */ |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */ |
} |
} |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
savm=oldm; |
savm=oldm; |
oldm=newm; |
oldm=newm; |
} /* end mult */ |
} /* end mult */ |
|
|
/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ |
/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ |
/* But now since version 0.9 we anticipate for bias at large stepm. |
/* But now since version 0.9 we anticipate for bias at large stepm. |
* If stepm is larger than one month (smallest stepm) and if the exact delay |
* If stepm is larger than one month (smallest stepm) and if the exact delay |
* (in months) between two waves is not a multiple of stepm, we rounded to |
* (in months) between two waves is not a multiple of stepm, we rounded to |
* the nearest (and in case of equal distance, to the lowest) interval but now |
* the nearest (and in case of equal distance, to the lowest) interval but now |
* we keep into memory the bias bh[mi][i] and also the previous matrix product |
* we keep into memory the bias bh[mi][i] and also the previous matrix product |
* (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the |
* (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the |
* probability in order to take into account the bias as a fraction of the way |
* probability in order to take into account the bias as a fraction of the way |
* from savm to out if bh is negative or even beyond if bh is positive. bh varies |
* from savm to out if bh is negative or even beyond if bh is positive. bh varies |
* -stepm/2 to stepm/2 . |
* -stepm/2 to stepm/2 . |
* For stepm=1 the results are the same as for previous versions of Imach. |
* For stepm=1 the results are the same as for previous versions of Imach. |
* For stepm > 1 the results are less biased than in previous versions. |
* For stepm > 1 the results are less biased than in previous versions. |
*/ |
*/ |
s1=s[mw[mi][i]][i]; |
s1=s[mw[mi][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
bbh=(double)bh[mi][i]/(double)stepm; |
bbh=(double)bh[mi][i]/(double)stepm; |
/* bias bh is positive if real duration |
/* bias bh is positive if real duration |
* is higher than the multiple of stepm and negative otherwise. |
* is higher than the multiple of stepm and negative otherwise. |
*/ |
*/ |
/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ |
/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ |
if( s2 > nlstate){ |
if( s2 > nlstate){ |
/* i.e. if s2 is a death state and if the date of death is known |
/* i.e. if s2 is a death state and if the date of death is known |
then the contribution to the likelihood is the probability to |
then the contribution to the likelihood is the probability to |
die between last step unit time and current step unit time, |
die between last step unit time and current step unit time, |
which is also equal to probability to die before dh |
which is also equal to probability to die before dh |
minus probability to die before dh-stepm . |
minus probability to die before dh-stepm . |
In version up to 0.92 likelihood was computed |
In version up to 0.92 likelihood was computed |
as if date of death was unknown. Death was treated as any other |
as if date of death was unknown. Death was treated as any other |
health state: the date of the interview describes the actual state |
health state: the date of the interview describes the actual state |
and not the date of a change in health state. The former idea was |
and not the date of a change in health state. The former idea was |
to consider that at each interview the state was recorded |
to consider that at each interview the state was recorded |
(healthy, disable or death) and IMaCh was corrected; but when we |
(healthy, disable or death) and IMaCh was corrected; but when we |
introduced the exact date of death then we should have modified |
introduced the exact date of death then we should have modified |
the contribution of an exact death to the likelihood. This new |
the contribution of an exact death to the likelihood. This new |
contribution is smaller and very dependent of the step unit |
contribution is smaller and very dependent of the step unit |
stepm. It is no more the probability to die between last interview |
stepm. It is no more the probability to die between last interview |
and month of death but the probability to survive from last |
and month of death but the probability to survive from last |
interview up to one month before death multiplied by the |
interview up to one month before death multiplied by the |
probability to die within a month. Thanks to Chris |
probability to die within a month. Thanks to Chris |
Jackson for correcting this bug. Former versions increased |
Jackson for correcting this bug. Former versions increased |
mortality artificially. The bad side is that we add another loop |
mortality artificially. The bad side is that we add another loop |
which slows down the processing. The difference can be up to 10% |
which slows down the processing. The difference can be up to 10% |
lower mortality. |
lower mortality. |
*/ |
*/ |
/* If, at the beginning of the maximization mostly, the |
/* If, at the beginning of the maximization mostly, the |
cumulative probability or probability to be dead is |
cumulative probability or probability to be dead is |
constant (ie = 1) over time d, the difference is equal to |
constant (ie = 1) over time d, the difference is equal to |
0. out[s1][3] = savm[s1][3]: probability, being at state |
0. out[s1][3] = savm[s1][3]: probability, being at state |
s1 at precedent wave, to be dead a month before current |
s1 at precedent wave, to be dead a month before current |
wave is equal to probability, being at state s1 at |
wave is equal to probability, being at state s1 at |
precedent wave, to be dead at mont of the current |
precedent wave, to be dead at mont of the current |
wave. Then the observed probability (that this person died) |
wave. Then the observed probability (that this person died) |
is null according to current estimated parameter. In fact, |
is null according to current estimated parameter. In fact, |
it should be very low but not zero otherwise the log go to |
it should be very low but not zero otherwise the log go to |
infinity. |
infinity. |
*/ |
*/ |
/* #ifdef INFINITYORIGINAL */ |
/* #ifdef INFINITYORIGINAL */ |
/* lli=log(out[s1][s2] - savm[s1][s2]); */ |
/* lli=log(out[s1][s2] - savm[s1][s2]); */ |
/* #else */ |
/* #else */ |
Line 3037 double func( double *x)
|
Line 3060 double func( double *x)
|
/* else */ |
/* else */ |
/* lli=log(out[s1][s2] - savm[s1][s2]); */ |
/* lli=log(out[s1][s2] - savm[s1][s2]); */ |
/* #endif */ |
/* #endif */ |
lli=log(out[s1][s2] - savm[s1][s2]); |
lli=log(out[s1][s2] - savm[s1][s2]); |
|
|
} else if ( s2==-1 ) { /* alive */ |
} else if ( s2==-1 ) { /* alive */ |
for (j=1,survp=0. ; j<=nlstate; j++) |
for (j=1,survp=0. ; j<=nlstate; j++) |
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
/*survp += out[s1][j]; */ |
/*survp += out[s1][j]; */ |
lli= log(survp); |
lli= log(survp); |
} |
} |
else if (s2==-4) { |
else if (s2==-4) { |
for (j=3,survp=0. ; j<=nlstate; j++) |
for (j=3,survp=0. ; j<=nlstate; j++) |
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
lli= log(survp); |
lli= log(survp); |
} |
} |
else if (s2==-5) { |
else if (s2==-5) { |
for (j=1,survp=0. ; j<=2; j++) |
for (j=1,survp=0. ; j<=2; j++) |
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; |
lli= log(survp); |
lli= log(survp); |
} |
} |
else{ |
else{ |
lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ |
lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ |
/* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ |
/* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ |
} |
} |
/*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ |
/*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ |
/*if(lli ==000.0)*/ |
/*if(lli ==000.0)*/ |
/*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ |
/*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ |
ipmx +=1; |
ipmx +=1; |
sw += weight[i]; |
sw += weight[i]; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
/* if (lli < log(mytinydouble)){ */ |
/* if (lli < log(mytinydouble)){ */ |
/* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */ |
/* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */ |
/* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */ |
/* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */ |
/* } */ |
/* } */ |
} /* end of wave */ |
} /* end of wave */ |
} /* end of individual */ |
} /* end of individual */ |
} else if(mle==2){ |
} else if(mle==2){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
} |
} |
for(d=0; d<=dh[mi][i]; d++){ |
for(d=0; d<=dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
cov[2]=agexact; |
cov[2]=agexact; |
if(nagesqr==1) |
if(nagesqr==1) |
cov[3]= agexact*agexact; |
cov[3]= agexact*agexact; |
for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
} |
} |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
savm=oldm; |
savm=oldm; |
oldm=newm; |
oldm=newm; |
} /* end mult */ |
} /* end mult */ |
|
|
s1=s[mw[mi][i]][i]; |
s1=s[mw[mi][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
bbh=(double)bh[mi][i]/(double)stepm; |
bbh=(double)bh[mi][i]/(double)stepm; |
lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ |
lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ |
ipmx +=1; |
ipmx +=1; |
sw += weight[i]; |
sw += weight[i]; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
} /* end of wave */ |
} /* end of wave */ |
} /* end of individual */ |
} /* end of individual */ |
} else if(mle==3){ /* exponential inter-extrapolation */ |
} else if(mle==3){ /* exponential inter-extrapolation */ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
} |
} |
for(d=0; d<dh[mi][i]; d++){ |
for(d=0; d<dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
cov[2]=agexact; |
cov[2]=agexact; |
if(nagesqr==1) |
if(nagesqr==1) |
cov[3]= agexact*agexact; |
cov[3]= agexact*agexact; |
for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
} |
} |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
savm=oldm; |
savm=oldm; |
oldm=newm; |
oldm=newm; |
} /* end mult */ |
} /* end mult */ |
|
|
s1=s[mw[mi][i]][i]; |
s1=s[mw[mi][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
bbh=(double)bh[mi][i]/(double)stepm; |
bbh=(double)bh[mi][i]/(double)stepm; |
lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ |
lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ |
ipmx +=1; |
ipmx +=1; |
sw += weight[i]; |
sw += weight[i]; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
} /* end of wave */ |
} /* end of wave */ |
} /* end of individual */ |
} /* end of individual */ |
}else if (mle==4){ /* ml=4 no inter-extrapolation */ |
}else if (mle==4){ /* ml=4 no inter-extrapolation */ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
} |
} |
for(d=0; d<dh[mi][i]; d++){ |
for(d=0; d<dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
cov[2]=agexact; |
cov[2]=agexact; |
if(nagesqr==1) |
if(nagesqr==1) |
cov[3]= agexact*agexact; |
cov[3]= agexact*agexact; |
for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
} |
} |
|
|
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
savm=oldm; |
savm=oldm; |
oldm=newm; |
oldm=newm; |
} /* end mult */ |
} /* end mult */ |
|
|
s1=s[mw[mi][i]][i]; |
s1=s[mw[mi][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
if( s2 > nlstate){ |
if( s2 > nlstate){ |
lli=log(out[s1][s2] - savm[s1][s2]); |
lli=log(out[s1][s2] - savm[s1][s2]); |
} else if ( s2==-1 ) { /* alive */ |
} else if ( s2==-1 ) { /* alive */ |
for (j=1,survp=0. ; j<=nlstate; j++) |
for (j=1,survp=0. ; j<=nlstate; j++) |
survp += out[s1][j]; |
survp += out[s1][j]; |
lli= log(survp); |
lli= log(survp); |
}else{ |
}else{ |
lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ |
lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ |
} |
} |
ipmx +=1; |
ipmx +=1; |
sw += weight[i]; |
sw += weight[i]; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
/* printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ |
/* printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */ |
} /* end of wave */ |
} /* end of wave */ |
} /* end of individual */ |
} /* end of individual */ |
}else{ /* ml=5 no inter-extrapolation no jackson =0.8a */ |
}else{ /* ml=5 no inter-extrapolation no jackson =0.8a */ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
} |
} |
for(d=0; d<dh[mi][i]; d++){ |
for(d=0; d<dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
cov[2]=agexact; |
cov[2]=agexact; |
if(nagesqr==1) |
if(nagesqr==1) |
cov[3]= agexact*agexact; |
cov[3]= agexact*agexact; |
for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
} |
} |
|
|
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
savm=oldm; |
savm=oldm; |
oldm=newm; |
oldm=newm; |
} /* end mult */ |
} /* end mult */ |
|
|
s1=s[mw[mi][i]][i]; |
s1=s[mw[mi][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
s2=s[mw[mi+1][i]][i]; |
lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ |
lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ |
ipmx +=1; |
ipmx +=1; |
sw += weight[i]; |
sw += weight[i]; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
/*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/ |
/*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/ |
} /* end of wave */ |
} /* end of wave */ |
} /* end of individual */ |
} /* end of individual */ |
} /* End of if */ |
} /* End of if */ |
for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; |
for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; |
/* printf("l1=%f l2=%f ",ll[1],ll[2]); */ |
/* printf("l1=%f l2=%f ",ll[1],ll[2]); */ |
l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ |
l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ |
return -l; |
return -l; |
} |
} |
|
|
/*************** log-likelihood *************/ |
/*************** log-likelihood *************/ |
Line 3248 double funcone( double *x)
|
Line 3271 double funcone( double *x)
|
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
ioffset=2+nagesqr+cptcovage; |
ioffset=2+nagesqr+cptcovage; |
/* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */ |
/* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */ |
for (k=1; k<=ncoveff;k++){ /* Simple and product covariates without age* products */ |
for (k=1; k<=ncoveff+nqfveff;k++){ /* Simple and product fixed covariates without age* products */ |
cov[++ioffset]=covar[Tvar[k]][i]; |
cov[++ioffset]=covar[Tvar[k]][i]; |
} |
} |
for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitatives Fixed covariates */ |
for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitative fixed covariates */ |
cov[++ioffset]=coqvar[iqv][i]; |
cov[++ioffset]=coqvar[Tvar[iqv]][i]; |
} |
} |
|
|
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ /* Varying with waves */ |
for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */ |
for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */ |
cov[ioffset+itv]=cotvar[mw[mi][i]][itv][i]; |
cov[ioffset+itv]=cotvar[mw[mi][i]][itv][i]; |
} |
} |
Line 3839 void pstamp(FILE *fichier)
|
Line 3862 void pstamp(FILE *fichier)
|
} |
} |
|
|
/************ Frequencies ********************/ |
/************ Frequencies ********************/ |
void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ |
void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ |
int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \ |
int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \ |
int firstpass, int lastpass, int stepm, int weightopt, char model[]) |
int firstpass, int lastpass, int stepm, int weightopt, char model[]) |
{ /* Some frequencies */ |
{ /* Some frequencies */ |
|
|
int i, m, jk, j1, bool, z1,j; |
int i, m, jk, j1, bool, z1,j; |
int iind=0, iage=0; |
int iind=0, iage=0; |
int mi; /* Effective wave */ |
int mi; /* Effective wave */ |
int first; |
int first; |
double ***freq; /* Frequencies */ |
double ***freq; /* Frequencies */ |
double *meanq; |
double *meanq; |
double **meanqt; |
double **meanqt; |
double *pp, **prop, *posprop, *pospropt; |
double *pp, **prop, *posprop, *pospropt; |
double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0; |
double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0; |
char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH]; |
char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH]; |
double agebegin, ageend; |
double agebegin, ageend; |
|
|
pp=vector(1,nlstate); |
pp=vector(1,nlstate); |
prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */ |
posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */ |
pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ |
pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ |
/* prop=matrix(1,nlstate,iagemin,iagemax+3); */ |
/* prop=matrix(1,nlstate,iagemin,iagemax+3); */ |
meanq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ |
meanq=vector(1,nqfveff); /* Number of Quantitative Fixed Variables Effective */ |
meanqt=matrix(1,lastpass,1,nqtveff); |
meanqt=matrix(1,lastpass,1,nqtveff); |
strcpy(fileresp,"P_"); |
strcpy(fileresp,"P_"); |
strcat(fileresp,fileresu); |
strcat(fileresp,fileresu); |
/*strcat(fileresphtm,fileresu);*/ |
/*strcat(fileresphtm,fileresu);*/ |
if((ficresp=fopen(fileresp,"w"))==NULL) { |
if((ficresp=fopen(fileresp,"w"))==NULL) { |
printf("Problem with prevalence resultfile: %s\n", fileresp); |
printf("Problem with prevalence resultfile: %s\n", fileresp); |
fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); |
fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); |
exit(0); |
exit(0); |
} |
} |
|
|
strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm")); |
strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm")); |
if((ficresphtm=fopen(fileresphtm,"w"))==NULL) { |
if((ficresphtm=fopen(fileresphtm,"w"))==NULL) { |
printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); |
printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); |
fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); |
fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno)); |
fflush(ficlog); |
fflush(ficlog); |
exit(70); |
exit(70); |
} |
} |
else{ |
else{ |
fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ |
fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
} |
} |
fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm); |
fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm); |
|
|
strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); |
strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); |
if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { |
if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { |
printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); |
printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); |
fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); |
fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno)); |
fflush(ficlog); |
fflush(ficlog); |
exit(70); |
exit(70); |
} |
} |
else{ |
else{ |
fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ |
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\ |
fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
} |
} |
fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr); |
fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr); |
|
|
freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE); |
j1=0; |
j1=0; |
|
|
j=ncoveff; |
j=ncoveff; |
if (cptcovn<1) {j=1;ncodemax[1]=1;} |
if (cptcovn<1) {j=1;ncodemax[1]=1;} |
|
|
first=1; |
first=1; |
|
|
/* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels: |
/* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels: |
reference=low_education V1=0,V2=0 |
reference=low_education V1=0,V2=0 |
med_educ V1=1 V2=0, |
med_educ V1=1 V2=0, |
high_educ V1=0 V2=1 |
high_educ V1=0 V2=1 |
Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff |
Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff |
*/ |
*/ |
|
|
for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination excluding varying and quantitatives */ |
for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination excluding varying and quantitatives */ |
posproptt=0.; |
posproptt=0.; |
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
scanf("%d", i);*/ |
scanf("%d", i);*/ |
for (i=-5; i<=nlstate+ndeath; i++) |
for (i=-5; i<=nlstate+ndeath; i++) |
for (jk=-5; jk<=nlstate+ndeath; jk++) |
for (jk=-5; jk<=nlstate+ndeath; jk++) |
for(m=iagemin; m <= iagemax+3; m++) |
for(m=iagemin; m <= iagemax+3; m++) |
freq[i][jk][m]=0; |
freq[i][jk][m]=0; |
|
|
for (i=1; i<=nlstate; i++) { |
for (i=1; i<=nlstate; i++) { |
for(m=iagemin; m <= iagemax+3; m++) |
for(m=iagemin; m <= iagemax+3; m++) |
prop[i][m]=0; |
prop[i][m]=0; |
posprop[i]=0; |
posprop[i]=0; |
pospropt[i]=0; |
pospropt[i]=0; |
} |
} |
for (z1=1; z1<= nqfveff; z1++) { |
for (z1=1; z1<= nqfveff; z1++) { |
meanq[z1]+=0.; |
meanq[z1]+=0.; |
for(m=1;m<=lastpass;m++){ |
for(m=1;m<=lastpass;m++){ |
meanqt[m][z1]=0.; |
meanqt[m][z1]=0.; |
} |
} |
} |
} |
|
|
dateintsum=0; |
|
k2cpt=0; |
|
/* For that comination of covariate j1, we count and print the frequencies */ |
|
for (iind=1; iind<=imx; iind++) { /* For each individual iind */ |
|
bool=1; |
|
if (nqfveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
|
for (z1=1; z1<= nqfveff; z1++) { |
|
meanq[z1]+=coqvar[Tvar[z1]][iind]; |
|
} |
|
for (z1=1; z1<=ncoveff; z1++) { |
|
/* if(Tvaraff[z1] ==-20){ */ |
|
/* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ |
|
/* }else if(Tvaraff[z1] ==-10){ */ |
|
/* /\* sumnew+=coqvar[z1][iind]; *\/ */ |
|
/* }else */ |
|
if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ |
|
/* Tests if this individual i responded to j1 (V4=1 V3=0) */ |
|
bool=0; |
|
/* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", |
|
bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), |
|
j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ |
|
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ |
|
} |
|
} /* end z1 */ |
|
} /* cptcovn > 0 */ |
|
|
|
if (bool==1){ /* We selected an individual iin satisfying combination j1 */ |
|
/* for(m=firstpass; m<=lastpass; m++){ */ |
|
for(mi=1; mi<wav[iind];mi++){ |
|
m=mw[mi][iind]; |
|
/* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind] |
|
and mw[mi+1][iind]. dh depends on stepm. */ |
|
agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/ |
|
ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */ |
|
if(m >=firstpass && m <=lastpass){ |
|
k2=anint[m][iind]+(mint[m][iind]/12.); |
|
/*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ |
|
if(agev[m][iind]==0) agev[m][iind]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ |
|
if(agev[m][iind]==1) agev[m][iind]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ |
|
if (s[m][iind]>0 && s[m][iind]<=nlstate) /* If status at wave m is known and a live state */ |
|
prop[s[m][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ |
|
if (m<lastpass) { |
|
/* if(s[m][iind]==4 && s[m+1][iind]==4) */ |
|
/* printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */ |
|
if(s[m][iind]==-1) |
|
printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.)); |
|
freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ |
|
/* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */ |
|
freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */ |
|
} |
|
} |
|
if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) { |
|
dateintsum=dateintsum+k2; |
|
k2cpt++; |
|
/* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ |
|
} |
|
/*}*/ |
|
} /* end m */ |
|
} /* end bool */ |
|
} /* end iind = 1 to imx */ |
|
/* prop[s][age] is feeded for any initial and valid live state as well as |
|
freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */ |
|
|
|
|
|
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
|
pstamp(ficresp); |
|
if (ncoveff>0) { |
|
fprintf(ficresp, "\n#********** Variable "); |
|
fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); |
|
fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); |
|
for (z1=1; z1<=ncoveff; z1++){ |
|
fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
|
fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
|
fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
|
} |
|
fprintf(ficresp, "**********\n#"); |
|
fprintf(ficresphtm, "**********</h3>\n"); |
|
fprintf(ficresphtmfr, "**********</h3>\n"); |
|
fprintf(ficlog, "\n#********** Variable "); |
|
for (z1=1; z1<=ncoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
|
fprintf(ficlog, "**********\n"); |
|
} |
|
fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">"); |
|
for(i=1; i<=nlstate;i++) { |
|
fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); |
|
fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i); |
|
} |
|
fprintf(ficresp, "\n"); |
|
fprintf(ficresphtm, "\n"); |
|
|
|
/* Header of frequency table by age */ |
|
fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">"); |
|
fprintf(ficresphtmfr,"<th>Age</th> "); |
|
for(jk=-1; jk <=nlstate+ndeath; jk++){ |
|
for(m=-1; m <=nlstate+ndeath; m++){ |
|
if(jk!=0 && m!=0) |
|
fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m); |
|
} |
|
} |
|
fprintf(ficresphtmfr, "\n"); |
|
|
|
/* For each age */ |
dateintsum=0; |
for(iage=iagemin; iage <= iagemax+3; iage++){ |
k2cpt=0; |
fprintf(ficresphtm,"<tr>"); |
/* For that comination of covariate j1, we count and print the frequencies */ |
if(iage==iagemax+1){ |
for (iind=1; iind<=imx; iind++) { /* For each individual iind */ |
fprintf(ficlog,"1"); |
bool=1; |
fprintf(ficresphtmfr,"<tr><th>0</th> "); |
if (nqfveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
}else if(iage==iagemax+2){ |
for (z1=1; z1<= nqfveff; z1++) { |
fprintf(ficlog,"0"); |
meanq[z1]+=coqvar[Tvar[z1]][iind]; |
fprintf(ficresphtmfr,"<tr><th>Unknown</th> "); |
} |
}else if(iage==iagemax+3){ |
for (z1=1; z1<=ncoveff; z1++) { |
fprintf(ficlog,"Total"); |
/* if(Tvaraff[z1] ==-20){ */ |
fprintf(ficresphtmfr,"<tr><th>Total</th> "); |
/* /\* sumnew+=cotvar[mw[mi][iind]][z1][iind]; *\/ */ |
}else{ |
/* }else if(Tvaraff[z1] ==-10){ */ |
if(first==1){ |
/* /\* sumnew+=coqvar[z1][iind]; *\/ */ |
first=0; |
/* }else */ |
printf("See log file for details...\n"); |
if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ |
} |
/* Tests if this individual i responded to j1 (V4=1 V3=0) */ |
fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage); |
bool=0; |
fprintf(ficlog,"Age %d", iage); |
/* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", |
} |
bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), |
for(jk=1; jk <=nlstate ; jk++){ |
j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ |
for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) |
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ |
pp[jk] += freq[jk][m][iage]; |
} |
} |
} /* end z1 */ |
for(jk=1; jk <=nlstate ; jk++){ |
} /* cptcovn > 0 */ |
for(m=-1, pos=0; m <=0 ; m++) |
|
pos += freq[jk][m][iage]; |
if (bool==1){ /* We selected an individual iin satisfying combination j1 */ |
if(pp[jk]>=1.e-10){ |
/* for(m=firstpass; m<=lastpass; m++){ */ |
if(first==1){ |
for(mi=1; mi<wav[iind];mi++){ |
printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
m=mw[mi][iind]; |
} |
/* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind] |
fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
and mw[mi+1][iind]. dh depends on stepm. */ |
}else{ |
agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/ |
if(first==1) |
ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */ |
printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); |
if(m >=firstpass && m <=lastpass){ |
fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); |
k2=anint[m][iind]+(mint[m][iind]/12.); |
} |
/*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/ |
} |
if(agev[m][iind]==0) agev[m][iind]=iagemax+1; /* All ages equal to 0 are in iagemax+1 */ |
|
if(agev[m][iind]==1) agev[m][iind]=iagemax+2; /* All ages equal to 1 are in iagemax+2 */ |
for(jk=1; jk <=nlstate ; jk++){ |
if (s[m][iind]>0 && s[m][iind]<=nlstate) /* If status at wave m is known and a live state */ |
/* posprop[jk]=0; */ |
prop[s[m][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ |
for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ |
if (m<lastpass) { |
pp[jk] += freq[jk][m][iage]; |
/* if(s[m][iind]==4 && s[m+1][iind]==4) */ |
} /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */ |
/* printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */ |
|
if(s[m][iind]==-1) |
for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){ |
printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.)); |
pos += pp[jk]; /* pos is the total number of transitions until this age */ |
freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ |
posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state |
/* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */ |
from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */ |
pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state |
} |
from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
} |
} |
if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) { |
for(jk=1; jk <=nlstate ; jk++){ |
dateintsum=dateintsum+k2; |
if(pos>=1.e-5){ |
k2cpt++; |
if(first==1) |
/* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ |
printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); |
} |
fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); |
/*}*/ |
}else{ |
} /* end m */ |
if(first==1) |
} /* end bool */ |
printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); |
} /* end iind = 1 to imx */ |
fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); |
/* prop[s][age] is feeded for any initial and valid live state as well as |
} |
freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */ |
if( iage <= iagemax){ |
|
if(pos>=1.e-5){ |
|
fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); |
/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/ |
fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); |
pstamp(ficresp); |
/*probs[iage][jk][j1]= pp[jk]/pos;*/ |
if (ncoveff>0) { |
/*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ |
fprintf(ficresp, "\n#********** Variable "); |
} |
fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); |
else{ |
fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); |
fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); |
for (z1=1; z1<=ncoveff; z1++){ |
fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta); |
fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
} |
fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
} |
fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
pospropt[jk] +=posprop[jk]; |
} |
} /* end loop jk */ |
fprintf(ficresp, "**********\n#"); |
/* pospropt=0.; */ |
fprintf(ficresphtm, "**********</h3>\n"); |
for(jk=-1; jk <=nlstate+ndeath; jk++){ |
fprintf(ficresphtmfr, "**********</h3>\n"); |
for(m=-1; m <=nlstate+ndeath; m++){ |
fprintf(ficlog, "\n#********** Variable "); |
if(freq[jk][m][iage] !=0 ) { /* minimizing output */ |
for (z1=1; z1<=ncoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
if(first==1){ |
fprintf(ficlog, "**********\n"); |
printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); |
} |
} |
fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">"); |
fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); |
for(i=1; i<=nlstate;i++) { |
} |
fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); |
if(jk!=0 && m!=0) |
fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i); |
fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]); |
} |
} |
fprintf(ficresp, "\n"); |
} /* end loop jk */ |
fprintf(ficresphtm, "\n"); |
posproptt=0.; |
|
for(jk=1; jk <=nlstate; jk++){ |
/* Header of frequency table by age */ |
posproptt += pospropt[jk]; |
fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">"); |
} |
fprintf(ficresphtmfr,"<th>Age</th> "); |
fprintf(ficresphtmfr,"</tr>\n "); |
for(jk=-1; jk <=nlstate+ndeath; jk++){ |
if(iage <= iagemax){ |
for(m=-1; m <=nlstate+ndeath; m++){ |
fprintf(ficresp,"\n"); |
if(jk!=0 && m!=0) |
fprintf(ficresphtm,"</tr>\n"); |
fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m); |
} |
} |
if(first==1) |
} |
printf("Others in log...\n"); |
fprintf(ficresphtmfr, "\n"); |
fprintf(ficlog,"\n"); |
|
} /* end loop age iage */ |
/* For each age */ |
fprintf(ficresphtm,"<tr><th>Tot</th>"); |
for(iage=iagemin; iage <= iagemax+3; iage++){ |
for(jk=1; jk <=nlstate ; jk++){ |
fprintf(ficresphtm,"<tr>"); |
if(posproptt < 1.e-5){ |
if(iage==iagemax+1){ |
fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt); |
fprintf(ficlog,"1"); |
}else{ |
fprintf(ficresphtmfr,"<tr><th>0</th> "); |
fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt); |
}else if(iage==iagemax+2){ |
} |
fprintf(ficlog,"0"); |
} |
fprintf(ficresphtmfr,"<tr><th>Unknown</th> "); |
fprintf(ficresphtm,"</tr>\n"); |
}else if(iage==iagemax+3){ |
fprintf(ficresphtm,"</table>\n"); |
fprintf(ficlog,"Total"); |
fprintf(ficresphtmfr,"</table>\n"); |
fprintf(ficresphtmfr,"<tr><th>Total</th> "); |
if(posproptt < 1.e-5){ |
}else{ |
fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1); |
if(first==1){ |
fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1); |
first=0; |
fprintf(ficres,"\n This combination (%d) is not valid and no result will be produced\n\n",j1); |
printf("See log file for details...\n"); |
invalidvarcomb[j1]=1; |
} |
}else{ |
fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage); |
fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1); |
fprintf(ficlog,"Age %d", iage); |
invalidvarcomb[j1]=0; |
} |
} |
for(jk=1; jk <=nlstate ; jk++){ |
fprintf(ficresphtmfr,"</table>\n"); |
for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++) |
} /* end selected combination of covariate j1 */ |
pp[jk] += freq[jk][m][iage]; |
dateintmean=dateintsum/k2cpt; |
} |
|
for(jk=1; jk <=nlstate ; jk++){ |
|
for(m=-1, pos=0; m <=0 ; m++) |
|
pos += freq[jk][m][iage]; |
|
if(pp[jk]>=1.e-10){ |
|
if(first==1){ |
|
printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
|
} |
|
fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); |
|
}else{ |
|
if(first==1) |
|
printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); |
|
fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk); |
|
} |
|
} |
|
|
|
for(jk=1; jk <=nlstate ; jk++){ |
|
/* posprop[jk]=0; */ |
|
for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */ |
|
pp[jk] += freq[jk][m][iage]; |
|
} /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */ |
|
|
|
for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){ |
|
pos += pp[jk]; /* pos is the total number of transitions until this age */ |
|
posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state |
|
from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
|
pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state |
|
from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */ |
|
} |
|
for(jk=1; jk <=nlstate ; jk++){ |
|
if(pos>=1.e-5){ |
|
if(first==1) |
|
printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); |
|
fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos); |
|
}else{ |
|
if(first==1) |
|
printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); |
|
fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk); |
|
} |
|
if( iage <= iagemax){ |
|
if(pos>=1.e-5){ |
|
fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); |
|
fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta); |
|
/*probs[iage][jk][j1]= pp[jk]/pos;*/ |
|
/*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/ |
|
} |
|
else{ |
|
fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta); |
|
fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta); |
|
} |
|
} |
|
pospropt[jk] +=posprop[jk]; |
|
} /* end loop jk */ |
|
/* pospropt=0.; */ |
|
for(jk=-1; jk <=nlstate+ndeath; jk++){ |
|
for(m=-1; m <=nlstate+ndeath; m++){ |
|
if(freq[jk][m][iage] !=0 ) { /* minimizing output */ |
|
if(first==1){ |
|
printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]); |
|
} |
|
fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]); |
|
} |
|
if(jk!=0 && m!=0) |
|
fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]); |
|
} |
|
} /* end loop jk */ |
|
posproptt=0.; |
|
for(jk=1; jk <=nlstate; jk++){ |
|
posproptt += pospropt[jk]; |
|
} |
|
fprintf(ficresphtmfr,"</tr>\n "); |
|
if(iage <= iagemax){ |
|
fprintf(ficresp,"\n"); |
|
fprintf(ficresphtm,"</tr>\n"); |
|
} |
|
if(first==1) |
|
printf("Others in log...\n"); |
|
fprintf(ficlog,"\n"); |
|
} /* end loop age iage */ |
|
fprintf(ficresphtm,"<tr><th>Tot</th>"); |
|
for(jk=1; jk <=nlstate ; jk++){ |
|
if(posproptt < 1.e-5){ |
|
fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt); |
|
}else{ |
|
fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt); |
|
} |
|
} |
|
fprintf(ficresphtm,"</tr>\n"); |
|
fprintf(ficresphtm,"</table>\n"); |
|
fprintf(ficresphtmfr,"</table>\n"); |
|
if(posproptt < 1.e-5){ |
|
fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1); |
|
fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1); |
|
fprintf(ficres,"\n This combination (%d) is not valid and no result will be produced\n\n",j1); |
|
invalidvarcomb[j1]=1; |
|
}else{ |
|
fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1); |
|
invalidvarcomb[j1]=0; |
|
} |
|
fprintf(ficresphtmfr,"</table>\n"); |
|
} /* end selected combination of covariate j1 */ |
|
dateintmean=dateintsum/k2cpt; |
|
|
fclose(ficresp); |
fclose(ficresp); |
fclose(ficresphtm); |
fclose(ficresphtm); |
fclose(ficresphtmfr); |
fclose(ficresphtmfr); |
free_vector(meanq,1,nqfveff); |
free_vector(meanq,1,nqfveff); |
free_matrix(meanqt,1,lastpass,1,nqtveff); |
free_matrix(meanqt,1,lastpass,1,nqtveff); |
free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE); |
free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE); |
free_vector(pospropt,1,nlstate); |
free_vector(pospropt,1,nlstate); |
free_vector(posprop,1,nlstate); |
free_vector(posprop,1,nlstate); |
free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE); |
free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE); |
free_vector(pp,1,nlstate); |
free_vector(pp,1,nlstate); |
/* End of freqsummary */ |
/* End of freqsummary */ |
} |
} |
|
|
/************ Prevalence ********************/ |
/************ Prevalence ********************/ |
void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) |
void prevalence(double ***probs, double agemin, double agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, int firstpass, int lastpass) |
Line 7590 int readdata(char datafile[], int firsto
|
Line 7613 int readdata(char datafile[], int firsto
|
return 1; |
return 1; |
} |
} |
coqvar[iv][i]=dval; |
coqvar[iv][i]=dval; |
|
covar[ncovcol+iv][i]=dval; /* including qvar in standard covar for performance reasons */ |
} |
} |
strcpy(line,stra); |
strcpy(line,stra); |
}/* end loop nqv */ |
}/* end loop nqv */ |
Line 7798 int decodemodel ( char model[], int last
|
Line 7822 int decodemodel ( char model[], int last
|
cptcovprod--; |
cptcovprod--; |
cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */ |
cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */ |
Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */ |
Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */ |
Typevar[k]=1; /* 2 for age product */ |
Typevar[k]=1; /* 1 for age product */ |
cptcovage++; /* Sums the number of covariates which include age as a product */ |
cptcovage++; /* Sums the number of covariates which include age as a product */ |
Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */ |
Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */ |
/*printf("stre=%s ", stre);*/ |
/*printf("stre=%s ", stre);*/ |
Line 7863 int decodemodel ( char model[], int last
|
Line 7887 int decodemodel ( char model[], int last
|
|
|
/* Decodemodel knows only the grammar (simple, product, age*) of the model but not what kind |
/* Decodemodel knows only the grammar (simple, product, age*) of the model but not what kind |
of variable (dummy vs quantitative, fixed vs time varying) is behind */ |
of variable (dummy vs quantitative, fixed vs time varying) is behind */ |
/* ncovcol= 1, nqv=1, ntv=2, nqtv= 1 = 5 possible variables data |
/* ncovcol= 1, nqv=1 | ntv=2, nqtv= 1 = 5 possible variables data: 2 fixed 3, varying |
model= V2 + V4 +V3 + V4*V3 + V5*age + V5 , V1 is not used saving its place |
model= V5 + V4 +V3 + V4*V3 + V5*age + V2 + V1*V2 + V1*age + V5*age, V1 is not used saving its place |
k = 1 2 3 4 5 6 |
k = 1 2 3 4 5 6 7 8 9 |
Tvar[k]= 2 4 3 1+1+2+1+1=6 5 5 |
Tvar[k]= 5 4 3 1+1+2+1+1=6 5 2 7 1 5 |
Typevar[k]=0 0 0 2 1 0 |
Typevar[k]= 0 0 0 2 1 0 2 1 1 |
|
Fixed[Tvar[k]]1 1 1 1 2 0 1 2 3 |
|
Dummy[Tvar[k]]1 0 0 0 2 1 1 2 3 |
*/ |
*/ |
/* Dispatching between quantitative and time varying covariates */ |
/* Dispatching between quantitative and time varying covariates */ |
|
/* If Tvar[k] >ncovcol it is a product */ |
/* Tvar[k] is the value n of Vn with n varying for 1 to nvcol, or p Vp=Vn*Vm for product */ |
/* Tvar[k] is the value n of Vn with n varying for 1 to nvcol, or p Vp=Vn*Vm for product */ |
|
/* Computing effective variables, ie used by the model, that is from the cptcovt variables */ |
for(k=1, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ |
for(k=1, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */ |
if (Tvar[k] <=ncovcol){ /* Simple fixed dummy covariatee */ |
if (Tvar[k] <=ncovcol && (Typevar[k]==0 || Typevar[k]==2)){ /* Simple or product fixed dummy covariatee */ |
|
Fixed[Tvar[k]]= 0; |
|
Dummy[Tvar[k]]= 0; |
ncoveff++; |
ncoveff++; |
}else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ |
}else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ |
|
Fixed[Tvar[k]]= 0; |
|
Dummy[Tvar[k]]= 1; |
nqfveff++; /* Only simple fixed quantitative variable */ |
nqfveff++; /* Only simple fixed quantitative variable */ |
}else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){ |
}else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 0; |
ntveff++; /* Only simple time varying dummy variable */ |
ntveff++; /* Only simple time varying dummy variable */ |
}else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv && Typevar[k]==0){ |
}else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){ |
nqtveff++;/* Only simple time varying quantitative variable */ |
if( Typevar[k]==0){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 1; |
|
nqtveff++;/* Only simple time varying quantitative variable */ |
|
} |
|
}else if (Typevar[k] == 2) { |
|
for(k1=1; k1 <= cptcovprodnoage; k1++){ |
|
if(Tvard[k1][1] <=ncovcol){ |
|
if(Tvard[k1][2] <=ncovcol){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 0; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv){ |
|
Fixed[Tvar[k]]= 0; |
|
Dummy[Tvar[k]]= 1; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 0; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 1; |
|
} |
|
}else if(Tvard[k1][1] <=ncovcol+nqv){ |
|
if(Tvard[k1][2] <=ncovcol){ |
|
Fixed[Tvar[k]]= 0; |
|
Dummy[Tvar[k]]= 1; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv){ |
|
Fixed[Tvar[k]]= 0; |
|
Dummy[Tvar[k]]= 1; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 1; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 1; |
|
} |
|
}else if(Tvard[k1][1] <=ncovcol+nqv+ntv){ |
|
if(Tvard[k1][2] <=ncovcol){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 1; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 1; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 0; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 1; |
|
} |
|
}else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){ |
|
if(Tvard[k1][2] <=ncovcol){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 1; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 1; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 1; |
|
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){ |
|
Fixed[Tvar[k]]= 1; |
|
Dummy[Tvar[k]]= 1; |
|
} |
|
}else{ |
|
printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); |
|
fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]); |
|
} |
|
} /* end k1 */ |
}else{ |
}else{ |
printf("Other types in effective covariates \n"); |
printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]); |
|
fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]); |
} |
} |
|
printf("Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[Tvar[k]],Dummy[Tvar[k]]); |
|
fprintf(ficlog,"Decodemodel, k=%d, Tvar[%d]=V%d,Typevar=%d, Fixed=%d, Dummy=%d\n",k, k,Tvar[k],Typevar[k],Fixed[Tvar[k]],Dummy[Tvar[k]]); |
} |
} |
|
printf("Model=%s\n\ |
|
Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ |
|
Fixed[Tvar[k]] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product \n\ |
|
Dummy[Tvar[k]] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); |
|
|
printf("ncoveff=%d, nqfveff=%d, ntveff=%d, nqtveff=%d, cptcovn=%d\n",ncoveff,nqfveff,ntveff,nqtveff,cptcovn); |
printf("ncoveff=%d, nqfveff=%d, ntveff=%d, nqtveff=%d, cptcovn=%d\n",ncoveff,nqfveff,ntveff,nqtveff,cptcovn); |
fprintf(ficlog,"ncoveff=%d, nqfveff=%d, ntveff=%d, nqtveff=%d, cptcovn=%d\n",ncoveff,nqfveff,ntveff,nqtveff,cptcovn); |
fprintf(ficlog,"ncoveff=%d, nqfveff=%d, ntveff=%d, nqtveff=%d, cptcovn=%d\n",ncoveff,nqfveff,ntveff,nqtveff,cptcovn); |
return (0); /* with covar[new additional covariate if product] and Tage if age */ |
return (0); /* with covar[new additional covariate if product] and Tage if age */ |
Line 9036 run imach with mle=-1 to get a correct t
|
Line 9144 run imach with mle=-1 to get a correct t
|
|
|
/* Scans npar lines */ |
/* Scans npar lines */ |
for(i=1; i <=npar; i++){ |
for(i=1; i <=npar; i++){ |
count=fscanf(ficpar,"%1d%1d%1d",&i1,&j1,&jk); |
count=fscanf(ficpar,"%1d%1d%d",&i1,&j1,&jk); |
if(count != 3){ |
if(count != 3){ |
printf("Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ |
printf("Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ |
This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ |
This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ |
Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); |
Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); |
fprintf(ficlog,"Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ |
fprintf(ficlog,"Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ |
This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ |
This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ |
Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); |
Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); |
exit(1); |
exit(1); |
}else{ |
}else{ |
if(mle==1) |
if(mle==1) |
printf("%1d%1d%1d",i1,j1,jk); |
printf("%1d%1d%d",i1,j1,jk); |
} |
} |
fprintf(ficlog,"%1d%1d%1d",i1,j1,jk); |
fprintf(ficlog,"%1d%1d%d",i1,j1,jk); |
fprintf(ficparo,"%1d%1d%1d",i1,j1,jk); |
fprintf(ficparo,"%1d%1d%d",i1,j1,jk); |
for(j=1; j <=i; j++){ |
for(j=1; j <=i; j++){ |
fscanf(ficpar," %le",&matcov[i][j]); |
fscanf(ficpar," %le",&matcov[i][j]); |
if(mle==1){ |
if(mle==1){ |
printf(" %.5le",matcov[i][j]); |
printf(" %.5le",matcov[i][j]); |
} |
} |
fprintf(ficlog," %.5le",matcov[i][j]); |
fprintf(ficlog," %.5le",matcov[i][j]); |
fprintf(ficparo," %.5le",matcov[i][j]); |
fprintf(ficparo," %.5le",matcov[i][j]); |
} |
} |
fscanf(ficpar,"\n"); |
fscanf(ficpar,"\n"); |
numlinepar++; |
numlinepar++; |
Line 9069 Please run with mle=-1 to get a correct
|
Line 9177 Please run with mle=-1 to get a correct
|
/* End of read covariance matrix npar lines */ |
/* End of read covariance matrix npar lines */ |
for(i=1; i <=npar; i++) |
for(i=1; i <=npar; i++) |
for(j=i+1;j<=npar;j++) |
for(j=i+1;j<=npar;j++) |
matcov[i][j]=matcov[j][i]; |
matcov[i][j]=matcov[j][i]; |
|
|
if(mle==1) |
if(mle==1) |
printf("\n"); |
printf("\n"); |
Line 9129 Please run with mle=-1 to get a correct
|
Line 9237 Please run with mle=-1 to get a correct
|
k=1 Tvar[1]=2 (from V2) |
k=1 Tvar[1]=2 (from V2) |
*/ |
*/ |
Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */ |
Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */ |
Typevar=ivector(1,NCOVMAX); /* -1 to 4 */ |
Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */ |
|
Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */ |
|
Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */ |
/* V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs). |
/* V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs). |
For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4, |
For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4, |
Tvar[4=age*V3] is 3 and 'age' is recorded in Tage. |
Tvar[4=age*V3] is 3 and 'age' is recorded in Tage. |
Line 9416 Interval (in months) between two waves:
|
Line 9526 Interval (in months) between two waves:
|
for (i=1; i<=imx; i++){ |
for (i=1; i<=imx; i++){ |
dcwave[i]=-1; |
dcwave[i]=-1; |
for (m=firstpass; m<=lastpass; m++) |
for (m=firstpass; m<=lastpass; m++) |
if (s[m][i]>nlstate) { |
if (s[m][i]>nlstate) { |
dcwave[i]=m; |
dcwave[i]=m; |
/* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/ |
/* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/ |
break; |
break; |
} |
} |
} |
} |
|
|
for (i=1; i<=imx; i++) { |
for (i=1; i<=imx; i++) { |
if (wav[i]>0){ |
if (wav[i]>0){ |
ageexmed[i]=agev[mw[1][i]][i]; |
ageexmed[i]=agev[mw[1][i]][i]; |
j=wav[i]; |
j=wav[i]; |
agecens[i]=1.; |
agecens[i]=1.; |
|
|
if (ageexmed[i]> 1 && wav[i] > 0){ |
if (ageexmed[i]> 1 && wav[i] > 0){ |
agecens[i]=agev[mw[j][i]][i]; |
agecens[i]=agev[mw[j][i]][i]; |
cens[i]= 1; |
cens[i]= 1; |
}else if (ageexmed[i]< 1) |
}else if (ageexmed[i]< 1) |
cens[i]= -1; |
cens[i]= -1; |
if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass) |
if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass) |
cens[i]=0 ; |
cens[i]=0 ; |
} |
} |
else cens[i]=-1; |
else cens[i]=-1; |
} |
} |
|
|
for (i=1;i<=NDIM;i++) { |
for (i=1;i<=NDIM;i++) { |
for (j=1;j<=NDIM;j++) |
for (j=1;j<=NDIM;j++) |
ximort[i][j]=(i == j ? 1.0 : 0.0); |
ximort[i][j]=(i == j ? 1.0 : 0.0); |
} |
} |
|
|
/*p[1]=0.0268; p[NDIM]=0.083;*/ |
/*p[1]=0.0268; p[NDIM]=0.083;*/ |
Line 10275 Please run with mle=-1 to get a correct
|
Line 10385 Please run with mle=-1 to get a correct
|
|
|
free_ivector(ncodemax,1,NCOVMAX); |
free_ivector(ncodemax,1,NCOVMAX); |
free_ivector(ncodemaxwundef,1,NCOVMAX); |
free_ivector(ncodemaxwundef,1,NCOVMAX); |
|
free_ivector(Dummy,-1,NCOVMAX); |
|
free_ivector(Fixed,-1,NCOVMAX); |
free_ivector(Typevar,-1,NCOVMAX); |
free_ivector(Typevar,-1,NCOVMAX); |
free_ivector(Tvar,1,NCOVMAX); |
free_ivector(Tvar,1,NCOVMAX); |
free_ivector(Tprod,1,NCOVMAX); |
free_ivector(Tprod,1,NCOVMAX); |