version 1.303, 2021/02/11 19:50:15
|
version 1.328, 2022/07/27 17:40:48
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
$State$ |
$State$ |
$Log$ |
$Log$ |
|
Revision 1.328 2022/07/27 17:40:48 brouard |
|
Summary: valgrind bug fixed by initializing to zero DummyV as well as Tage |
|
|
|
Revision 1.327 2022/07/27 14:47:35 brouard |
|
Summary: Still a problem for one-step probabilities in case of quantitative variables |
|
|
|
Revision 1.326 2022/07/26 17:33:55 brouard |
|
Summary: some test with nres=1 |
|
|
|
Revision 1.325 2022/07/25 14:27:23 brouard |
|
Summary: r30 |
|
|
|
* imach.c (Module): Error cptcovn instead of nsd in bmij (was |
|
coredumped, revealed by Feiuno, thank you. |
|
|
|
Revision 1.324 2022/07/23 17:44:26 brouard |
|
*** empty log message *** |
|
|
|
Revision 1.323 2022/07/22 12:30:08 brouard |
|
* imach.c (Module): Output of Wald test in the htm file and not only in the log. |
|
|
|
Revision 1.322 2022/07/22 12:27:48 brouard |
|
* imach.c (Module): Output of Wald test in the htm file and not only in the log. |
|
|
|
Revision 1.321 2022/07/22 12:04:24 brouard |
|
Summary: r28 |
|
|
|
* imach.c (Module): Output of Wald test in the htm file and not only in the log. |
|
|
|
Revision 1.320 2022/06/02 05:10:11 brouard |
|
*** empty log message *** |
|
|
|
Revision 1.319 2022/06/02 04:45:11 brouard |
|
* imach.c (Module): Adding the Wald tests from the log to the main |
|
htm for better display of the maximum likelihood estimators. |
|
|
|
Revision 1.318 2022/05/24 08:10:59 brouard |
|
* imach.c (Module): Some attempts to find a bug of wrong estimates |
|
of confidencce intervals with product in the equation modelC |
|
|
|
Revision 1.317 2022/05/15 15:06:23 brouard |
|
* imach.c (Module): Some minor improvements |
|
|
|
Revision 1.316 2022/05/11 15:11:31 brouard |
|
Summary: r27 |
|
|
|
Revision 1.315 2022/05/11 15:06:32 brouard |
|
*** empty log message *** |
|
|
|
Revision 1.314 2022/04/13 17:43:09 brouard |
|
* imach.c (Module): Adding link to text data files |
|
|
|
Revision 1.313 2022/04/11 15:57:42 brouard |
|
* imach.c (Module): Error in rewriting the 'r' file with yearsfproj or yearsbproj fixed |
|
|
|
Revision 1.312 2022/04/05 21:24:39 brouard |
|
*** empty log message *** |
|
|
|
Revision 1.311 2022/04/05 21:03:51 brouard |
|
Summary: Fixed quantitative covariates |
|
|
|
Fixed covariates (dummy or quantitative) |
|
with missing values have never been allowed but are ERRORS and |
|
program quits. Standard deviations of fixed covariates were |
|
wrongly computed. Mean and standard deviations of time varying |
|
covariates are still not computed. |
|
|
|
Revision 1.310 2022/03/17 08:45:53 brouard |
|
Summary: 99r25 |
|
|
|
Improving detection of errors: result lines should be compatible with |
|
the model. |
|
|
|
Revision 1.309 2021/05/20 12:39:14 brouard |
|
Summary: Version 0.99r24 |
|
|
|
Revision 1.308 2021/03/31 13:11:57 brouard |
|
Summary: Version 0.99r23 |
|
|
|
|
|
* imach.c (Module): Still bugs in the result loop. Thank to Holly Benett |
|
|
|
Revision 1.307 2021/03/08 18:11:32 brouard |
|
Summary: 0.99r22 fixed bug on result: |
|
|
|
Revision 1.306 2021/02/20 15:44:02 brouard |
|
Summary: Version 0.99r21 |
|
|
|
* imach.c (Module): Fix bug on quitting after result lines! |
|
(Module): Version 0.99r21 |
|
|
|
Revision 1.305 2021/02/20 15:28:30 brouard |
|
* imach.c (Module): Fix bug on quitting after result lines! |
|
|
|
Revision 1.304 2021/02/12 11:34:20 brouard |
|
* imach.c (Module): The use of a Windows BOM (huge) file is now an error |
|
|
Revision 1.303 2021/02/11 19:50:15 brouard |
Revision 1.303 2021/02/11 19:50:15 brouard |
* (Module): imach.c Someone entered 'results:' instead of 'result:'. Now it is an error which is printed. |
* (Module): imach.c Someone entered 'results:' instead of 'result:'. Now it is an error which is printed. |
|
|
Line 790
|
Line 887
|
|
|
The same imach parameter file can be used but the option for mle should be -3. |
The same imach parameter file can be used but the option for mle should be -3. |
|
|
Agnès, who wrote this part of the code, tried to keep most of the |
Agnès, who wrote this part of the code, tried to keep most of the |
former routines in order to include the new code within the former code. |
former routines in order to include the new code within the former code. |
|
|
The output is very simple: only an estimate of the intercept and of |
The output is very simple: only an estimate of the intercept and of |
Line 969 Important routines
|
Line 1066 Important routines
|
- Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities) |
- 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. |
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 |
- 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 |
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. |
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 |
from the European Union. |
from the European Union. |
It is copyrighted identically to a GNU software product, ie programme and |
It is copyrighted identically to a GNU software product, ie programme and |
Line 1039 Important routines
|
Line 1136 Important routines
|
#define POWELLNOF3INFF1TEST /* Skip test */ |
#define POWELLNOF3INFF1TEST /* Skip test */ |
/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ |
/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ |
/* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */ |
/* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */ |
|
/* #define FLATSUP *//* Suppresses directions where likelihood is flat */ |
|
|
#include <math.h> |
#include <math.h> |
#include <stdio.h> |
#include <stdio.h> |
Line 1105 typedef struct {
|
Line 1203 typedef struct {
|
#define NINTERVMAX 8 |
#define NINTERVMAX 8 |
#define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ |
#define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ |
#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ |
#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ |
#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ |
#define NCOVMAX 30 /**< Maximum number of covariates used in the model, including generated covariates V1*V2 or V1*age */ |
#define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
#define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
/*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/ |
/*#define decodtabm(h,k,cptcoveff)= (h <= (1<<cptcoveff)?(((h-1) >> (k-1)) & 1) +1 : -1)*/ |
#define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 |
#define decodtabm(h,k,cptcoveff) (((h-1) >> (k-1)) & 1) +1 |
Line 1133 typedef struct {
|
Line 1231 typedef struct {
|
/* $State$ */ |
/* $State$ */ |
#include "version.h" |
#include "version.h" |
char version[]=__IMACH_VERSION__; |
char version[]=__IMACH_VERSION__; |
char copyright[]="May 2019,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020"; |
char copyright[]="July 2022,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022"; |
char fullversion[]="$Revision$ $Date$"; |
char fullversion[]="$Revision$ $Date$"; |
char strstart[80]; |
char strstart[80]; |
char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
Line 1315 double ***cotvar; /* Time varying covari
|
Line 1413 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 */ |
/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
/* Some documentation */ |
/*k 1 2 3 4 5 6 7 8 9 */ |
/* Design original data |
/*Tvar[k]= 5 4 3 6 5 2 7 1 1 */ |
* V1 V2 V3 V4 V5 V6 V7 V8 Weight ddb ddth d1st s1 V9 V10 V11 V12 s2 V9 V10 V11 V12 |
/* Tndvar[k] 1 2 3 4 5 */ |
* < ncovcol=6 > nqv=2 (V7 V8) dv dv dv qtv dv dv dvv qtv |
/*TDvar 4 3 6 7 1 */ /* For outputs only; combination of dummies fixed or varying */ |
* ntv=3 nqtv=1 |
/* Tns[k] 1 2 2 4 5 */ /* Number of single cova */ |
* cptcovn number of covariates (not including constant and age) = # of + plus 1 = 10+1=11 |
/* TvarsD[k] 1 2 3 */ /* Number of single dummy cova */ |
* For time varying covariate, quanti or dummies |
/* TvarsDind 2 3 9 */ /* position K of single dummy cova */ |
* cotqvar[wav][iv(1 to nqtv)][i]= [1][12][i]=(V12) quanti |
/* TvarsQ[k] 1 2 */ /* Number of single quantitative cova */ |
* cotvar[wav][ntv+iv][i]= [3+(1 to nqtv)][i]=(V12) quanti |
/* TvarsQind 1 6 */ /* position K of single quantitative cova */ |
* cotvar[wav][iv(1 to ntv)][i]= [1][1][i]=(V9) dummies at wav 1 |
/* Tprod[i]=k 4 7 */ |
* cotvar[wav][iv(1 to ntv)][i]= [1][2][i]=(V10) dummies at wav 1 |
/* Tage[i]=k 5 8 */ |
* covar[k,i], value of kth fixed covariate dummy or quanti : |
/* */ |
* covar[1][i]= (V1), covar[4][i]=(V4), covar[8][i]=(V8) |
|
* Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 + V9 + V9*age + V10 |
|
* k= 1 2 3 4 5 6 7 8 9 10 11 |
|
*/ |
|
/* According to the model, more columns can be added to covar by the product of covariates */ |
|
/* ncovcol=1(Males=0 Females=1) nqv=1(raedyrs) ntv=2(withoutiadl=0 withiadl=1, witoutadl=0 withoutadl=1) nqtv=1(bmi) nlstate=3 ndeath=1 |
|
# States 1=Coresidence, 2 Living alone, 3 Institution |
|
# V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi |
|
*/ |
|
/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
|
/* k 1 2 3 4 5 6 7 8 9 */ |
|
/*Typevar[k]= 0 0 0 2 1 0 2 1 0 *//*0 for simple covariate (dummy, quantitative,*/ |
|
/* fixed or varying), 1 for age product, 2 for*/ |
|
/* product */ |
|
/*Dummy[k]= 1 0 0 1 3 1 1 2 0 *//*Dummy[k] 0=dummy (0 1), 1 quantitative */ |
|
/*(single or product without age), 2 dummy*/ |
|
/* with age product, 3 quant with age product*/ |
|
/*Tvar[k]= 5 4 3 6 5 2 7 1 1 */ |
|
/* nsd 1 2 3 */ /* Counting single dummies covar fixed or tv */ |
|
/*TvarsD[nsd] 4 3 1 */ /* ID of single dummy cova fixed or timevary*/ |
|
/*TvarsDind[k] 2 3 9 */ /* position K of single dummy cova */ |
|
/* nsq 1 2 */ /* Counting single quantit tv */ |
|
/* TvarsQ[k] 5 2 */ /* Number of single quantitative cova */ |
|
/* TvarsQind 1 6 */ /* position K of single quantitative cova */ |
|
/* Tprod[i]=k 1 2 */ /* Position in model of the ith prod without age */ |
|
/* cptcovage 1 2 */ /* Counting cov*age in the model equation */ |
|
/* Tage[cptcovage]=k 5 8 */ /* Position in the model of ith cov*age */ |
|
/* Tvard[1][1]@4={4,3,1,2} V4*V3 V1*V2 */ /* Position in model of the ith prod without age */ |
|
/* TvarF TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 ID of fixed covariates or product V2, V1*V2, V1 */ |
|
/* TvarFind; TvarFind[1]=6, TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod) */ |
/* Type */ |
/* Type */ |
/* V 1 2 3 4 5 */ |
/* V 1 2 3 4 5 */ |
/* F F V V V */ |
/* F F V V V */ |
Line 1338 int *TvarsDind;
|
Line 1465 int *TvarsDind;
|
int *TvarsQ; |
int *TvarsQ; |
int *TvarsQind; |
int *TvarsQind; |
|
|
#define MAXRESULTLINES 10 |
#define MAXRESULTLINESPONE 10+1 |
int nresult=0; |
int nresult=0; |
int parameterline=0; /* # of the parameter (type) line */ |
int parameterline=0; /* # of the parameter (type) line */ |
int TKresult[MAXRESULTLINES]; |
int TKresult[MAXRESULTLINESPONE]; |
int Tresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */ |
int Tresult[MAXRESULTLINESPONE][NCOVMAX];/* For dummy variable , value (output) */ |
int Tinvresult[MAXRESULTLINES][NCOVMAX];/* For dummy variable , value (output) */ |
int Tinvresult[MAXRESULTLINESPONE][NCOVMAX];/* For dummy variable , value (output) */ |
int Tvresult[MAXRESULTLINES][NCOVMAX]; /* For dummy variable , variable # (output) */ |
int Tvresult[MAXRESULTLINESPONE][NCOVMAX]; /* For dummy variable , variable # (output) */ |
double Tqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */ |
double Tqresult[MAXRESULTLINESPONE][NCOVMAX]; /* For quantitative variable , value (output) */ |
double Tqinvresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , value (output) */ |
double Tqinvresult[MAXRESULTLINESPONE][NCOVMAX]; /* For quantitative variable , value (output) */ |
int Tvqresult[MAXRESULTLINES][NCOVMAX]; /* For quantitative variable , variable # (output) */ |
int Tvqresult[MAXRESULTLINESPONE][NCOVMAX]; /* For quantitative variable , variable # (output) */ |
|
|
|
/* ncovcol=1(Males=0 Females=1) nqv=1(raedyrs) ntv=2(withoutiadl=0 withiadl=1, witoutadl=0 withoutadl=1) nqtv=1(bmi) nlstate=3 ndeath=1 |
|
# States 1=Coresidence, 2 Living alone, 3 Institution |
|
# V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi |
|
*/ |
/* int *TDvar; /\**< TDvar[1]=4, TDvarF[2]=3, TDvar[3]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */ |
/* int *TDvar; /\**< TDvar[1]=4, TDvarF[2]=3, TDvar[3]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 *\/ */ |
int *TvarF; /**< TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
int *TvarF; /**< TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
int *TvarFind; /**< TvarFind[1]=6, TvarFind[2]=7, Tvarind[3]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
int *TvarFind; /**< TvarFind[1]=6, TvarFind[2]=7, Tvarind[3]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
Line 1550 char *cutl(char *blocc, char *alocc, cha
|
Line 1681 char *cutl(char *blocc, char *alocc, cha
|
{ |
{ |
/* cuts string in into blocc and alocc where blocc ends before FIRST occurence of char 'occ' |
/* cuts string in into blocc and alocc where blocc ends before FIRST occurence of char 'occ' |
and alocc starts after first occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2') |
and alocc starts after first occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2') |
gives blocc="abcdef" and alocc="ghi2j". |
gives alocc="abcdef" and blocc="ghi2j". |
If occ is not found blocc is null and alocc is equal to in. Returns blocc |
If occ is not found blocc is null and alocc is equal to in. Returns blocc |
*/ |
*/ |
char *s, *t; |
char *s, *t; |
Line 1832 char *subdirf(char fileres[])
|
Line 1963 char *subdirf(char fileres[])
|
/*************** function subdirf2 ***********/ |
/*************** function subdirf2 ***********/ |
char *subdirf2(char fileres[], char *preop) |
char *subdirf2(char fileres[], char *preop) |
{ |
{ |
|
/* Example subdirf2(optionfilefiname,"FB_") with optionfilefiname="texte", result="texte/FB_texte" |
|
Errors in subdirf, 2, 3 while printing tmpout is |
|
rewritten within the same printf. Workaround: many printfs */ |
/* Caution optionfilefiname is hidden */ |
/* Caution optionfilefiname is hidden */ |
strcpy(tmpout,optionfilefiname); |
strcpy(tmpout,optionfilefiname); |
strcat(tmpout,"/"); |
strcat(tmpout,"/"); |
Line 2203 void linmin(double p[], double xi[], int
|
Line 2336 void linmin(double p[], double xi[], int
|
#endif |
#endif |
#ifdef LINMINORIGINAL |
#ifdef LINMINORIGINAL |
#else |
#else |
if(fb == fx){ /* Flat function in the direction */ |
if(fb == fx){ /* Flat function in the direction */ |
xmin=xx; |
xmin=xx; |
*flat=1; |
*flat=1; |
}else{ |
}else{ |
*flat=0; |
*flat=0; |
#endif |
#endif |
/*Flat mnbrak2 shift (*ax=0.000000000000, *fa=51626.272983130431), (*bx=-1.618034000000, *fb=51590.149499362531), (*cx=-4.236068025156, *fc=51590.149499362531) */ |
/*Flat mnbrak2 shift (*ax=0.000000000000, *fa=51626.272983130431), (*bx=-1.618034000000, *fb=51590.149499362531), (*cx=-4.236068025156, *fc=51590.149499362531) */ |
Line 2264 void linmin(double p[], double xi[], int
|
Line 2397 void linmin(double p[], double xi[], int
|
|
|
/*************** powell ************************/ |
/*************** powell ************************/ |
/* |
/* |
Minimization of a function func of n variables. Input consists of an initial starting point |
Minimization of a function func of n variables. Input consists in an initial starting point |
p[1..n] ; an initial matrix xi[1..n][1..n] , whose columns contain the initial set of di- |
p[1..n] ; an initial matrix xi[1..n][1..n] whose columns contain the initial set of di- |
rections (usually the n unit vectors); and ftol , the fractional tolerance in the function value |
rections (usually the n unit vectors); and ftol, the fractional tolerance in the function value |
such that failure to decrease by more than this amount on one iteration signals doneness. On |
such that failure to decrease by more than this amount in one iteration signals doneness. On |
output, p is set to the best point found, xi is the then-current direction set, fret is the returned |
output, p is set to the best point found, xi is the then-current direction set, fret is the returned |
function value at p , and iter is the number of iterations taken. The routine linmin is used. |
function value at p , and iter is the number of iterations taken. The routine linmin is used. |
*/ |
*/ |
Line 2292 void powell(double p[], double **xi, int
|
Line 2425 void powell(double p[], double **xi, int
|
double fp,fptt; |
double fp,fptt; |
double *xits; |
double *xits; |
int niterf, itmp; |
int niterf, itmp; |
#ifdef LINMINORIGINAL |
|
#else |
|
|
|
flatdir=ivector(1,n); |
|
for (j=1;j<=n;j++) flatdir[j]=0; |
|
#endif |
|
|
|
pt=vector(1,n); |
pt=vector(1,n); |
ptt=vector(1,n); |
ptt=vector(1,n); |
Line 2307 void powell(double p[], double **xi, int
|
Line 2434 void powell(double p[], double **xi, int
|
for (j=1;j<=n;j++) pt[j]=p[j]; |
for (j=1;j<=n;j++) pt[j]=p[j]; |
rcurr_time = time(NULL); |
rcurr_time = time(NULL); |
for (*iter=1;;++(*iter)) { |
for (*iter=1;;++(*iter)) { |
fp=(*fret); /* From former iteration or initial value */ |
|
ibig=0; |
ibig=0; |
del=0.0; |
del=0.0; |
rlast_time=rcurr_time; |
rlast_time=rcurr_time; |
/* (void) gettimeofday(&curr_time,&tzp); */ |
/* (void) gettimeofday(&curr_time,&tzp); */ |
rcurr_time = time(NULL); |
rcurr_time = time(NULL); |
curr_time = *localtime(&rcurr_time); |
curr_time = *localtime(&rcurr_time); |
printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); |
printf("\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); |
fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); |
fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); |
/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */ |
/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */ |
|
fp=(*fret); /* From former iteration or initial value */ |
for (i=1;i<=n;i++) { |
for (i=1;i<=n;i++) { |
fprintf(ficrespow," %.12lf", p[i]); |
fprintf(ficrespow," %.12lf", p[i]); |
} |
} |
Line 2421 void powell(double p[], double **xi, int
|
Line 2548 void powell(double p[], double **xi, int
|
/* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */ |
/* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */ |
/* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */ |
/* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */ |
/* New value of last point Pn is not computed, P(n-1) */ |
/* New value of last point Pn is not computed, P(n-1) */ |
for(j=1;j<=n;j++) { |
for(j=1;j<=n;j++) { |
if(flatdir[j] >0){ |
if(flatdir[j] >0){ |
printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); |
printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); |
fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); |
fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); |
} |
|
/* printf("\n"); */ |
|
/* fprintf(ficlog,"\n"); */ |
|
} |
} |
|
/* printf("\n"); */ |
|
/* fprintf(ficlog,"\n"); */ |
|
} |
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */ |
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */ |
if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */ |
if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */ |
/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ |
/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ |
Line 2466 void powell(double p[], double **xi, int
|
Line 2593 void powell(double p[], double **xi, int
|
} |
} |
#endif |
#endif |
|
|
#ifdef LINMINORIGINAL |
|
#else |
|
free_ivector(flatdir,1,n); |
|
#endif |
|
free_vector(xit,1,n); |
free_vector(xit,1,n); |
free_vector(xits,1,n); |
free_vector(xits,1,n); |
free_vector(ptt,1,n); |
free_vector(ptt,1,n); |
Line 2583 void powell(double p[], double **xi, int
|
Line 2706 void powell(double p[], double **xi, int
|
} |
} |
printf("\n"); |
printf("\n"); |
fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
|
#ifdef FLATSUP |
|
free_vector(xit,1,n); |
|
free_vector(xits,1,n); |
|
free_vector(ptt,1,n); |
|
free_vector(pt,1,n); |
|
return; |
|
#endif |
} |
} |
#endif |
#endif |
printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); |
printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); |
Line 2667 void powell(double p[], double **xi, int
|
Line 2797 void powell(double p[], double **xi, int
|
newm=savm; |
newm=savm; |
/* Covariates have to be included here again */ |
/* Covariates have to be included here again */ |
cov[2]=agefin; |
cov[2]=agefin; |
if(nagesqr==1) |
if(nagesqr==1){ |
cov[3]= agefin*agefin;; |
cov[3]= agefin*agefin; |
|
} |
for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ |
for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ |
/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ |
/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ |
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; |
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; |
|
/* cov[++k1]=nbcode[TvarsD[k]][codtabm(ij,k)]; */ |
/* printf("prevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ |
/* printf("prevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ |
} |
} |
for (k=1; k<=nsq;k++) { /* For single varying covariates only */ |
for (k=1; k<=nsq;k++) { /* For single varying covariates only */ |
/* Here comes the value of quantitative after renumbering k with single quantitative covariates */ |
/* Here comes the value of quantitative after renumbering k with single quantitative covariates */ |
cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; |
cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; |
|
/* cov[++k1]=Tqresult[nres][k]; */ |
/* printf("prevalim Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */ |
/* printf("prevalim Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */ |
} |
} |
for (k=1; k<=cptcovage;k++){ /* For product with age */ |
for (k=1; k<=cptcovage;k++){ /* For product with age */ |
if(Dummy[Tvar[Tage[k]]]){ |
if(Dummy[Tage[k]]==2){ /* dummy with age */ |
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
} else{ |
/* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ |
cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; |
} else if(Dummy[Tage[k]]==3){ /* quantitative with age */ |
|
cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; |
|
/* cov[++k1]=Tqresult[nres][k]; */ |
} |
} |
/* printf("prevalim Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ |
/* printf("prevalim Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ |
} |
} |
Line 2692 void powell(double p[], double **xi, int
|
Line 2827 void powell(double p[], double **xi, int
|
if(Dummy[Tvard[k][1]==0]){ |
if(Dummy[Tvard[k][1]==0]){ |
if(Dummy[Tvard[k][2]==0]){ |
if(Dummy[Tvard[k][2]==0]){ |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; |
|
/* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */ |
}else{ |
}else{ |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; |
|
/* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */ |
} |
} |
}else{ |
}else{ |
if(Dummy[Tvard[k][2]==0]){ |
if(Dummy[Tvard[k][2]==0]){ |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; |
|
/* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; */ |
}else{ |
}else{ |
cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; |
cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; |
|
/* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; */ |
} |
} |
} |
} |
} |
} |
Line 2708 void powell(double p[], double **xi, int
|
Line 2847 void powell(double p[], double **xi, int
|
/*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/ |
/*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/ |
/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ |
/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ |
/* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */ |
/* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */ |
/* age and covariate values of ij are in 'cov' */ |
/* age and covariate values of ij are in 'cov' */ |
out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* Bug Valgrind */ |
out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* Bug Valgrind */ |
|
|
savm=oldm; |
savm=oldm; |
Line 2748 void powell(double p[], double **xi, int
|
Line 2887 void powell(double p[], double **xi, int
|
if(!first){ |
if(!first){ |
first=1; |
first=1; |
printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d). Others in log file only...\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); |
printf("Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d). Others in log file only...\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); |
|
fprintf(ficlog, "Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d).\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); |
|
}else if (first >=1 && first <10){ |
|
fprintf(ficlog, "Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d).\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); |
|
first++; |
|
}else if (first ==10){ |
|
fprintf(ficlog, "Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d).\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); |
|
printf("Warning: the stable prevalence dit not converge. This warning came too often, IMaCh will stop notifying, even in its log file. Look at the graphs to appreciate the non convergence.\n"); |
|
fprintf(ficlog,"Warning: the stable prevalence no convergence; too many cases, giving up noticing, even in log file\n"); |
|
first++; |
} |
} |
fprintf(ficlog, "Warning: the stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.d years and %d loops. Try to lower 'ftolpl'. Youngest age to start was %d=(%d-%d).\n", (int)age, maxmax, ftolpl, *ncvyear, ncvloop, (int)(agefin+stepm/YEARM), (int)(age-stepm/YEARM), (int)delaymax); |
|
|
|
/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ |
/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */ |
free_vector(min,1,nlstate); |
free_vector(min,1,nlstate); |
Line 2823 void powell(double p[], double **xi, int
|
Line 2970 void powell(double p[], double **xi, int
|
/* newm points to the allocated table savm passed by the function it can be written, savm could be reallocated */ |
/* newm points to the allocated table savm passed by the function it can be written, savm could be reallocated */ |
/* Covariates have to be included here again */ |
/* Covariates have to be included here again */ |
cov[2]=agefin; |
cov[2]=agefin; |
if(nagesqr==1) |
if(nagesqr==1){ |
cov[3]= agefin*agefin;; |
cov[3]= agefin*agefin;; |
|
} |
for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ |
for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ |
/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ |
/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ |
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; |
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; |
Line 2845 void powell(double p[], double **xi, int
|
Line 2993 void powell(double p[], double **xi, int
|
/* /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; *\/ */ |
/* /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; *\/ */ |
/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */ |
/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */ |
for (k=1; k<=cptcovage;k++){ /* For product with age */ |
for (k=1; k<=cptcovage;k++){ /* For product with age */ |
if(Dummy[Tvar[Tage[k]]]){ |
/* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age *\/ ERROR ???*/ |
|
if(Dummy[Tage[k]]== 2){ /* dummy with age */ |
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
} else{ |
} else if(Dummy[Tage[k]]== 3){ /* quantitative with age */ |
cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; |
cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; |
} |
} |
/* printf("prevalim Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ |
/* printf("prevalim Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ |
} |
} |
Line 2918 void powell(double p[], double **xi, int
|
Line 3067 void powell(double p[], double **xi, int
|
|
|
maxmax=0.; |
maxmax=0.; |
for(i=1; i<=nlstate; i++){ |
for(i=1; i<=nlstate; i++){ |
meandiff[i]=(max[i]-min[i])/(max[i]+min[i])*2.; /* mean difference for each column */ |
meandiff[i]=(max[i]-min[i])/(max[i]+min[i])*2.; /* mean difference for each column, could be nan! */ |
maxmax=FMAX(maxmax,meandiff[i]); |
maxmax=FMAX(maxmax,meandiff[i]); |
/* printf("Back age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, i, meandiff[i],(int)agefin, i, max[i], i, min[i],maxmax); */ |
/* printf("Back age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, i, meandiff[i],(int)agefin, i, max[i], i, min[i],maxmax); */ |
} /* i loop */ |
} /* i loop */ |
Line 3004 double **pmij(double **ps, double *cov,
|
Line 3153 double **pmij(double **ps, double *cov,
|
ps[i][i]=1./(s1+1.); |
ps[i][i]=1./(s1+1.); |
/* Computing other pijs */ |
/* Computing other pijs */ |
for(j=1; j<i; j++) |
for(j=1; j<i; j++) |
ps[i][j]= exp(ps[i][j])*ps[i][i]; |
ps[i][j]= exp(ps[i][j])*ps[i][i];/* Bug valgrind */ |
for(j=i+1; j<=nlstate+ndeath; j++) |
for(j=i+1; j<=nlstate+ndeath; j++) |
ps[i][j]= exp(ps[i][j])*ps[i][i]; |
ps[i][j]= exp(ps[i][j])*ps[i][i]; |
/* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */ |
/* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */ |
Line 3052 double **pmij(double **ps, double *cov,
|
Line 3201 double **pmij(double **ps, double *cov,
|
doldm=ddoldms; /* global pointers */ |
doldm=ddoldms; /* global pointers */ |
dnewm=ddnewms; |
dnewm=ddnewms; |
dsavm=ddsavms; |
dsavm=ddsavms; |
|
|
|
/* Debug */ |
|
/* printf("Bmij ij=%d, cov[2}=%f\n", ij, cov[2]); */ |
agefin=cov[2]; |
agefin=cov[2]; |
/* Bx = Diag(w_x) P_x Diag(Sum_i w^i_x p^ij_x */ |
/* Bx = Diag(w_x) P_x Diag(Sum_i w^i_x p^ij_x */ |
/* bmij *//* age is cov[2], ij is included in cov, but we need for |
/* bmij *//* age is cov[2], ij is included in cov, but we need for |
Line 3060 double **pmij(double **ps, double *cov,
|
Line 3211 double **pmij(double **ps, double *cov,
|
/* dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ |
/* dsavm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ |
|
|
/* P_x */ |
/* P_x */ |
pmmij=pmij(pmmij,cov,ncovmodel,x,nlstate); /*This is forward probability from agefin to agefin + stepm */ |
pmmij=pmij(pmmij,cov,ncovmodel,x,nlstate); /*This is forward probability from agefin to agefin + stepm *//* Bug valgrind */ |
/* outputs pmmij which is a stochastic matrix in row */ |
/* outputs pmmij which is a stochastic matrix in row */ |
|
|
/* Diag(w_x) */ |
/* Diag(w_x) */ |
Line 3272 double ***hpxij(double ***po, int nhstep
|
Line 3423 double ***hpxij(double ***po, int nhstep
|
cov[1]=1.; |
cov[1]=1.; |
agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */ |
agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */ |
cov[2]=agexact; |
cov[2]=agexact; |
if(nagesqr==1) |
if(nagesqr==1){ |
cov[3]= agexact*agexact; |
cov[3]= agexact*agexact; |
|
} |
for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ |
for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ |
/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ |
/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */ |
|
/* codtabm(ij,k) (1 & (ij-1) >> (k-1))+1 */ |
|
/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
|
/* k 1 2 3 4 5 6 7 8 9 */ |
|
/*Tvar[k]= 5 4 3 6 5 2 7 1 1 */ |
|
/* nsd 1 2 3 */ /* Counting single dummies covar fixed or tv */ |
|
/*TvarsD[nsd] 4 3 1 */ /* ID of single dummy cova fixed or timevary*/ |
|
/*TvarsDind[k] 2 3 9 */ /* position K of single dummy cova */ |
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; |
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; |
/* printf("hpxij Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ |
/* printf("hpxij Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ |
} |
} |
for (k=1; k<=nsq;k++) { /* For single varying covariates only */ |
for (k=1; k<=nsq;k++) { /* For single varying covariates only */ |
/* Here comes the value of quantitative after renumbering k with single quantitative covariates */ |
/* Here comes the value of quantitative after renumbering k with single quantitative covariates */ |
cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; |
cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; |
/* printf("hPxij Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */ |
/* printf("hPxij Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */ |
} |
} |
for (k=1; k<=cptcovage;k++){ |
for (k=1; k<=cptcovage;k++){ /* For product with age V1+V1*age +V4 +age*V3 */ |
if(Dummy[Tvar[Tage[k]]]){ |
/* 1+2 Tage[1]=2 TVar[2]=1 Dummy[2]=2, Tage[2]=4 TVar[4]=3 Dummy[4]=3 quant*/ |
|
/* */ |
|
if(Dummy[Tage[k]]== 2){ /* dummy with age */ |
|
/* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age *\/ */ |
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
} else{ |
} else if(Dummy[Tage[k]]== 3){ /* quantitative with age */ |
cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; |
cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; |
} |
} |
/* printf("hPxij Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ |
/* printf("hPxij Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ |
} |
} |
for (k=1; k<=cptcovprod;k++){ /* */ |
for (k=1; k<=cptcovprod;k++){ /* For product without age */ |
/* printf("hPxij Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */ |
/* printf("hPxij Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */ |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; |
/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */ |
|
if(Dummy[Tvard[k][1]==0]){ |
|
if(Dummy[Tvard[k][2]==0]){ |
|
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; |
|
}else{ |
|
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; |
|
} |
|
}else{ |
|
if(Dummy[Tvard[k][2]==0]){ |
|
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; |
|
}else{ |
|
cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; |
|
} |
|
} |
} |
} |
/* for (k=1; k<=cptcovn;k++) */ |
/* for (k=1; k<=cptcovn;k++) */ |
/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ |
/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ |
Line 3306 double ***hpxij(double ***po, int nhstep
|
Line 3481 double ***hpxij(double ***po, int nhstep
|
|
|
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ |
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ |
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ |
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ |
/* right multiplication of oldm by the current matrix */ |
/* right multiplication of oldm by the current matrix */ |
out=matprod2(newm,oldm,1,nlstate+ndeath,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)); |
pmij(pmmij,cov,ncovmodel,x,nlstate)); |
/* if((int)age == 70){ */ |
/* if((int)age == 70){ */ |
Line 3376 double ***hbxij(double ***po, int nhstep
|
Line 3551 double ***hbxij(double ***po, int nhstep
|
cov[1]=1.; |
cov[1]=1.; |
agexact=age-( (h-1)*hstepm + (d) )*stepm/YEARM; /* age just before transition, d or d-1? */ |
agexact=age-( (h-1)*hstepm + (d) )*stepm/YEARM; /* age just before transition, d or d-1? */ |
/* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ |
/* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */ |
|
/* Debug */ |
|
/* printf("hBxij age=%lf, agexact=%lf\n", age, agexact); */ |
cov[2]=agexact; |
cov[2]=agexact; |
if(nagesqr==1) |
if(nagesqr==1) |
cov[3]= agexact*agexact; |
cov[3]= agexact*agexact; |
for (k=1; k<=cptcovn;k++){ |
for (k=1; k<=nsd;k++){ /* For single dummy covariates only *//* cptcovn error */ |
/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ |
/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */ |
/* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */ |
/* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */ |
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)]; |
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];/* Bug valgrind */ |
/* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ |
/* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */ |
} |
} |
for (k=1; k<=nsq;k++) { /* For single varying covariates only */ |
for (k=1; k<=nsq;k++) { /* For single varying covariates only */ |
Line 3390 double ***hbxij(double ***po, int nhstep
|
Line 3567 double ***hbxij(double ***po, int nhstep
|
cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; |
cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; |
/* printf("hPxij Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */ |
/* printf("hPxij Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */ |
} |
} |
for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 */ |
for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 *//* For product with age */ |
if(Dummy[Tvar[Tage[k]]]){ |
/* if(Dummy[Tvar[Tage[k]]]== 2){ /\* dummy with age error!!!*\/ */ |
|
if(Dummy[Tage[k]]== 2){ /* dummy with age */ |
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
} else{ |
} else if(Dummy[Tage[k]]== 3){ /* quantitative with age */ |
cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; |
cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; |
} |
} |
/* printf("hBxij Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ |
/* printf("hBxij Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */ |
} |
} |
for (k=1; k<=cptcovprod;k++){ /* Useless because included in cptcovn */ |
for (k=1; k<=cptcovprod;k++){ /* Useless because included in cptcovn */ |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; |
|
if(Dummy[Tvard[k][1]==0]){ |
|
if(Dummy[Tvard[k][2]==0]){ |
|
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; |
|
}else{ |
|
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; |
|
} |
|
}else{ |
|
if(Dummy[Tvard[k][2]==0]){ |
|
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; |
|
}else{ |
|
cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; |
|
} |
|
} |
} |
} |
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ |
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ |
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ |
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/ |
Line 3409 double ***hbxij(double ***po, int nhstep
|
Line 3600 double ***hbxij(double ***po, int nhstep
|
/* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */ |
/* out=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); */ |
out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\ |
out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent,ij),\ |
1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); |
1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);/* Bug valgrind */ |
/* if((int)age == 70){ */ |
/* if((int)age == 70){ */ |
/* printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */ |
/* printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */ |
/* for(i=1; i<=nlstate+ndeath; i++) { */ |
/* for(i=1; i<=nlstate+ndeath; i++) { */ |
Line 3495 double func( double *x)
|
Line 3686 double func( double *x)
|
*/ |
*/ |
ioffset=2+nagesqr ; |
ioffset=2+nagesqr ; |
/* Fixed */ |
/* Fixed */ |
for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */ |
for (k=1; k<=ncovf;k++){ /* For each fixed covariate dummu or quant or prod */ |
cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/ |
/* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */ |
|
/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
|
/* TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 ID of fixed covariates or product V2, V1*V2, V1 */ |
|
/* TvarFind; TvarFind[1]=6, TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod) */ |
|
cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (TvarFind[1]=6)*/ |
|
/* V1*V2 (7) TvarFind[2]=7, TvarFind[3]=9 */ |
} |
} |
/* 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 5, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]=6 |
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 |
Line 3511 double func( double *x)
|
Line 3707 double func( double *x)
|
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(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/ |
for(k=1; k <= ncovv ; k++){ /* Varying covariates in the model (single and product but no age )"V5+V4+V3+V4*V3+V5*age+V1*age+V1" +TvarVind 1,2,3,4(V4*V3) Tvar[1]@7{5, 4, 3, 6, 5, 1, 1 ; 6 because the created covar is after V5 and is 6, minus 1+1, 3,2,1,4 positions in cotvar*/ |
/* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */ |
/* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; but where is the crossproduct? */ |
cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; |
cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; |
} |
} |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (ii=1;ii<=nlstate+ndeath;ii++) |
Line 3527 double func( double *x)
|
Line 3723 double func( double *x)
|
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++) { |
if(!FixedV[Tvar[Tage[kk]]]) |
if(!FixedV[Tvar[Tage[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 */ |
else |
else |
cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact; |
cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][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)); |
Line 3638 double func( double *x)
|
Line 3834 double func( double *x)
|
} /* 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]; |
ioffset=2+nagesqr ; |
|
for (k=1; k<=ncovf;k++) |
|
cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i]; |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
|
for(k=1; k <= ncovv ; k++){ |
|
cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; |
|
} |
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); |
Line 3817 double funcone( double *x)
|
Line 4018 double funcone( double *x)
|
/* Fixed */ |
/* Fixed */ |
/* 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 fixed Dummy covariates without age* products *\/ */ |
/* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */ |
for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */ |
for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */ |
cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/ |
cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/ |
/* cov[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i]; */ |
/* cov[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i]; */ |
/* cov[2+6]=covar[Tvar[6]][i]; */ |
/* cov[2+6]=covar[Tvar[6]][i]; */ |
Line 3996 void likelione(FILE *ficres,double p[],
|
Line 4197 void likelione(FILE *ficres,double p[],
|
|
|
void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double [])) |
void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double [])) |
{ |
{ |
int i,j, iter=0; |
int i,j,k, jk, jkk=0, iter=0; |
double **xi; |
double **xi; |
double fret; |
double fret; |
double fretone; /* Only one call to likelihood */ |
double fretone; /* Only one call to likelihood */ |
Line 4030 void mlikeli(FILE *ficres,double p[], in
|
Line 4231 void mlikeli(FILE *ficres,double p[], in
|
if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); |
if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); |
fprintf(ficrespow,"\n"); |
fprintf(ficrespow,"\n"); |
#ifdef POWELL |
#ifdef POWELL |
|
#ifdef LINMINORIGINAL |
|
#else /* LINMINORIGINAL */ |
|
|
|
flatdir=ivector(1,npar); |
|
for (j=1;j<=npar;j++) flatdir[j]=0; |
|
#endif /*LINMINORIGINAL */ |
|
|
|
#ifdef FLATSUP |
|
powell(p,xi,npar,ftol,&iter,&fret,flatdir,func); |
|
/* reorganizing p by suppressing flat directions */ |
|
for(i=1, jk=1; i <=nlstate; i++){ |
|
for(k=1; k <=(nlstate+ndeath); k++){ |
|
if (k != i) { |
|
printf("%d%d flatdir[%d]=%d",i,k,jk, flatdir[jk]); |
|
if(flatdir[jk]==1){ |
|
printf(" To be skipped %d%d flatdir[%d]=%d ",i,k,jk, flatdir[jk]); |
|
} |
|
for(j=1; j <=ncovmodel; j++){ |
|
printf("%12.7f ",p[jk]); |
|
jk++; |
|
} |
|
printf("\n"); |
|
} |
|
} |
|
} |
|
/* skipping */ |
|
/* for(i=1, jk=1, jkk=1;(flatdir[jk]==0)&& (i <=nlstate); i++){ */ |
|
for(i=1, jk=1, jkk=1;i <=nlstate; i++){ |
|
for(k=1; k <=(nlstate+ndeath); k++){ |
|
if (k != i) { |
|
printf("%d%d flatdir[%d]=%d",i,k,jk, flatdir[jk]); |
|
if(flatdir[jk]==1){ |
|
printf(" To be skipped %d%d flatdir[%d]=%d jk=%d p[%d] ",i,k,jk, flatdir[jk],jk, jk); |
|
for(j=1; j <=ncovmodel; jk++,j++){ |
|
printf(" p[%d]=%12.7f",jk, p[jk]); |
|
/*q[jjk]=p[jk];*/ |
|
} |
|
}else{ |
|
printf(" To be kept %d%d flatdir[%d]=%d jk=%d q[%d]=p[%d] ",i,k,jk, flatdir[jk],jk, jkk, jk); |
|
for(j=1; j <=ncovmodel; jk++,jkk++,j++){ |
|
printf(" p[%d]=%12.7f=q[%d]",jk, p[jk],jkk); |
|
/*q[jjk]=p[jk];*/ |
|
} |
|
} |
|
printf("\n"); |
|
} |
|
fflush(stdout); |
|
} |
|
} |
|
powell(p,xi,npar,ftol,&iter,&fret,flatdir,func); |
|
#else /* FLATSUP */ |
powell(p,xi,npar,ftol,&iter,&fret,func); |
powell(p,xi,npar,ftol,&iter,&fret,func); |
#endif |
#endif /* FLATSUP */ |
|
|
|
#ifdef LINMINORIGINAL |
|
#else |
|
free_ivector(flatdir,1,npar); |
|
#endif /* LINMINORIGINAL*/ |
|
#endif /* POWELL */ |
|
|
#ifdef NLOPT |
#ifdef NLOPT |
#ifdef NEWUOA |
#ifdef NEWUOA |
Line 4059 void mlikeli(FILE *ficres,double p[], in
|
Line 4317 void mlikeli(FILE *ficres,double p[], in
|
} |
} |
nlopt_destroy(opt); |
nlopt_destroy(opt); |
#endif |
#endif |
|
#ifdef FLATSUP |
|
/* npared = npar -flatd/ncovmodel; */ |
|
/* xired= matrix(1,npared,1,npared); */ |
|
/* paramred= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); */ |
|
/* powell(pred,xired,npared,ftol,&iter,&fret,flatdir,func); */ |
|
/* free_matrix(xire,1,npared,1,npared); */ |
|
#else /* FLATSUP */ |
|
#endif /* FLATSUP */ |
free_matrix(xi,1,npar,1,npar); |
free_matrix(xi,1,npar,1,npar); |
fclose(ficrespow); |
fclose(ficrespow); |
printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); |
printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); |
Line 4510 void freqsummary(char fileres[], double
|
Line 4776 void freqsummary(char fileres[], double
|
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 and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm); |
fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies (weight=%d) and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm, weightopt); |
|
|
strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); |
strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm")); |
if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { |
if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) { |
Line 4520 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4786 Title=%s <br>Datafile=%s Firstpass=%d La
|
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 of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr); |
fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>(weight=%d) frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr,weightopt); |
|
|
y= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE); |
y= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE); |
x= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE); |
x= vector(iagemin-AGEMARGE,iagemax+4+AGEMARGE); |
Line 4672 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4938 Title=%s <br>Datafile=%s Firstpass=%d La
|
if(s[m][iind]==-1) |
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.)); |
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)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */ |
for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean */ |
for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean on known values only */ |
idq[z1]=idq[z1]+weight[iind]; |
if(!isnan(covar[ncovcol+z1][iind])){ |
meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /* Computes mean of quantitative with selected filter */ |
idq[z1]=idq[z1]+weight[iind]; |
stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /* *weight[iind];*/ /* Computes mean of quantitative with selected filter */ |
meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /* Computes mean of quantitative with selected filter */ |
|
/* stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; *//*error*/ |
|
stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]; /* *weight[iind];*/ /* Computes mean of quantitative with selected filter */ |
|
} |
} |
} |
/* if((int)agev[m][iind] == 55) */ |
/* if((int)agev[m][iind] == 55) */ |
/* printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */ |
/* printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */ |
Line 4699 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 4968 Title=%s <br>Datafile=%s Firstpass=%d La
|
/* } */ |
/* } */ |
} /* end bool */ |
} /* end bool */ |
} /* end iind = 1 to imx */ |
} /* end iind = 1 to imx */ |
/* prop[s][age] is feeded for any initial and valid live state as well as |
/* prop[s][age] is fed 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 */ |
freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */ |
|
|
|
|
Line 4738 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 5007 Title=%s <br>Datafile=%s Firstpass=%d La
|
Printing means of quantitative variables if any |
Printing means of quantitative variables if any |
*/ |
*/ |
for (z1=1; z1<= nqfveff; z1++) { |
for (z1=1; z1<= nqfveff; z1++) { |
fprintf(ficlog,"Mean of fixed quantitative variable V%d on %.0f individuals sum=%f", ncovcol+z1, idq[z1], meanq[z1]); |
fprintf(ficlog,"Mean of fixed quantitative variable V%d on %.3g (weighted) individuals sum=%f", ncovcol+z1, idq[z1], meanq[z1]); |
fprintf(ficlog,", mean=%.3g\n",meanq[z1]/idq[z1]); |
fprintf(ficlog,", mean=%.3g\n",meanq[z1]/idq[z1]); |
if(weightopt==1){ |
if(weightopt==1){ |
printf(" Weighted mean and standard deviation of"); |
printf(" Weighted mean and standard deviation of"); |
fprintf(ficlog," Weighted mean and standard deviation of"); |
fprintf(ficlog," Weighted mean and standard deviation of"); |
fprintf(ficresphtmfr," Weighted mean and standard deviation of"); |
fprintf(ficresphtmfr," Weighted mean and standard deviation of"); |
} |
} |
printf(" fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); |
/* mu = \frac{w x}{\sum w} |
fprintf(ficlog," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); |
var = \frac{\sum w (x-mu)^2}{\sum w} = \frac{w x^2}{\sum w} - mu^2 |
fprintf(ficresphtmfr," fixed quantitative variable V%d on %.0f representatives of the population : %6.3g (%6.3g)<p>\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt((stdq[z1]-meanq[z1]*meanq[z1]/idq[z1])/idq[z1])); |
*/ |
|
printf(" fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1])); |
|
fprintf(ficlog," fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1])); |
|
fprintf(ficresphtmfr," fixed quantitative variable V%d on %.3g (weighted) representatives of the population : %8.5g (%8.5g)<p>\n", ncovcol+z1, idq[z1],meanq[z1]/idq[z1], sqrt(stdq[z1]/idq[z1]-meanq[z1]*meanq[z1]/idq[z1]/idq[z1])); |
} |
} |
/* for (z1=1; z1<= nqtveff; z1++) { */ |
/* for (z1=1; z1<= nqtveff; z1++) { */ |
/* for(m=1;m<=lastpass;m++){ */ |
/* for(m=1;m<=lastpass;m++){ */ |
Line 5250 void concatwav(int wav[], int **dh, int
|
Line 5522 void concatwav(int wav[], int **dh, int
|
for(i=1; i<=imx; i++){ /* For simple cases and if state is death */ |
for(i=1; i<=imx; i++){ /* For simple cases and if state is death */ |
mi=0; /* First valid wave */ |
mi=0; /* First valid wave */ |
mli=0; /* Last valid wave */ |
mli=0; /* Last valid wave */ |
m=firstpass; |
m=firstpass; /* Loop on waves */ |
while(s[m][i] <= nlstate){ /* a live state */ |
while(s[m][i] <= nlstate){ /* a live state or unknown state */ |
if(m >firstpass && s[m][i]==s[m-1][i] && mint[m][i]==mint[m-1][i] && anint[m][i]==anint[m-1][i]){/* Two succesive identical information on wave m */ |
if(m >firstpass && s[m][i]==s[m-1][i] && mint[m][i]==mint[m-1][i] && anint[m][i]==anint[m-1][i]){/* Two succesive identical information on wave m */ |
mli=m-1;/* mw[++mi][i]=m-1; */ |
mli=m-1;/* mw[++mi][i]=m-1; */ |
}else if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */ |
}else if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */ |
mw[++mi][i]=m; |
mw[++mi][i]=m; /* Valid wave: incrementing mi and updating mi; mw[mi] is the wave number of mi_th valid transition */ |
mli=m; |
mli=m; |
} /* else might be a useless wave -1 and mi is not incremented and mw[mi] not updated */ |
} /* else might be a useless wave -1 and mi is not incremented and mw[mi] not updated */ |
if(m < lastpass){ /* m < lastpass, standard case */ |
if(m < lastpass){ /* m < lastpass, standard case */ |
m++; /* mi gives the "effective" current wave, m the current wave, go to next wave by incrementing m */ |
m++; /* mi gives the "effective" current wave, m the current wave, go to next wave by incrementing m */ |
} |
} |
else{ /* m >= lastpass, eventual special issue with warning */ |
else{ /* m = lastpass, eventual special issue with warning */ |
#ifdef UNKNOWNSTATUSNOTCONTRIBUTING |
#ifdef UNKNOWNSTATUSNOTCONTRIBUTING |
break; |
break; |
#else |
#else |
if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ |
if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){ /* no death date and known date of interview, case -2 (vital status unknown is warned later */ |
if(firsthree == 0){ |
if(firsthree == 0){ |
printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p_{%d%d} .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath); |
printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p_{%d%d} .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath); |
firsthree=1; |
firsthree=1; |
|
}else if(firsthree >=1 && firsthree < 10){ |
|
fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p_{%d%d} .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath); |
|
firsthree++; |
|
}else if(firsthree == 10){ |
|
printf("Information, too many Information flags: no more reported to log either\n"); |
|
fprintf(ficlog,"Information, too many Information flags: no more reported to log either\n"); |
|
firsthree++; |
|
}else{ |
|
firsthree++; |
} |
} |
fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as 1-p_{%d%d} .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m, s[m][i], nlstate+ndeath); |
mw[++mi][i]=m; /* Valid transition with unknown status */ |
mw[++mi][i]=m; |
|
mli=m; |
mli=m; |
} |
} |
if(s[m][i]==-2){ /* Vital status is really unknown */ |
if(s[m][i]==-2){ /* Vital status is really unknown */ |
nbwarn++; |
nbwarn++; |
if((int)anint[m][i] == 9999){ /* Has the vital status really been verified? */ |
if((int)anint[m][i] == 9999){ /* Has the vital status really been verified?not a transition */ |
printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); |
printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); |
fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); |
fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m); |
} |
} |
Line 5301 void concatwav(int wav[], int **dh, int
|
Line 5581 void concatwav(int wav[], int **dh, int
|
#ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE |
#ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE |
else if ((int) andc[i] != 9999) { /* Date of death is known */ |
else if ((int) andc[i] != 9999) { /* Date of death is known */ |
if ((int)anint[m][i]!= 9999) { /* date of last interview is known */ |
if ((int)anint[m][i]!= 9999) { /* date of last interview is known */ |
if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */ |
if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* month of death occured before last wave month and status should have been death instead of -1 */ |
nbwarn++; |
nbwarn++; |
if(firstfiv==0){ |
if(firstfiv==0){ |
printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); |
printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d, interviewed on %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); |
firstfiv=1; |
firstfiv=1; |
}else{ |
}else{ |
fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); |
fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d, interviewed on %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); |
} |
} |
}else{ /* Death occured afer last wave potential bias */ |
s[m][i]=nlstate+1; /* Fixing the status as death. Be careful if multiple death states */ |
|
}else{ /* Month of Death occured afer last wave month, potential bias */ |
nberr++; |
nberr++; |
if(firstwo==0){ |
if(firstwo==0){ |
printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d with status %d. Potential bias if other individuals are still alive on this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictitious wave at the date of last vital status scan, with a dead status. See documentation\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); |
firstwo=1; |
firstwo=1; |
} |
} |
fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictive wave at the date of last vital status scan, with a dead status or alive but unknown state status (-1). See documentation\n\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d with status %d. Potential bias if other individuals are still alive on this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood. Please add a new fictitious wave at the date of last vital status scan, with a dead status. See documentation\n\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); |
} |
} |
}else{ /* if date of interview is unknown */ |
}else{ /* if date of interview is unknown */ |
/* death is known but not confirmed by death status at any wave */ |
/* death is known but not confirmed by death status at any wave */ |
if(firstfour==0){ |
if(firstfour==0){ |
printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
printf("Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d with status %d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); |
firstfour=1; |
firstfour=1; |
} |
} |
fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m ); |
fprintf(ficlog,"Error! Death for individual %ld line=%d occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d with status %d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m ); |
} |
} |
} /* end if date of death is known */ |
} /* end if date of death is known */ |
#endif |
#endif |
wav[i]=mi; /* mi should be the last effective wave (or mli) */ |
wav[i]=mi; /* mi should be the last effective wave (or mli), */ |
/* wav[i]=mw[mi][i]; */ |
/* wav[i]=mw[mi][i]; */ |
if(mi==0){ |
if(mi==0){ |
nbwarn++; |
nbwarn++; |
if(first==0){ |
if(first==0){ |
Line 5342 void concatwav(int wav[], int **dh, int
|
Line 5623 void concatwav(int wav[], int **dh, int
|
} /* End individuals */ |
} /* End individuals */ |
/* wav and mw are no more changed */ |
/* wav and mw are no more changed */ |
|
|
|
printf("Information, you have to check %d informations which haven't been logged!\n",firsthree); |
|
fprintf(ficlog,"Information, you have to check %d informations which haven't been logged!\n",firsthree); |
|
|
|
|
for(i=1; i<=imx; i++){ |
for(i=1; i<=imx; i++){ |
for(mi=1; mi<wav[i];mi++){ |
for(mi=1; mi<wav[i];mi++){ |
if (stepm <=0) |
if (stepm <=0) |
Line 5463 void concatwav(int wav[], int **dh, int
|
Line 5747 void concatwav(int wav[], int **dh, int
|
if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ |
if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */ |
switch(Fixed[k]) { |
switch(Fixed[k]) { |
case 0: /* Testing on fixed dummy covariate, simple or product of fixed */ |
case 0: /* Testing on fixed dummy covariate, simple or product of fixed */ |
|
modmaxcovj=0; |
|
modmincovj=0; |
for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ |
for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ |
ij=(int)(covar[Tvar[k]][i]); |
ij=(int)(covar[Tvar[k]][i]); |
/* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i |
/* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i |
Line 5476 void concatwav(int wav[], int **dh, int
|
Line 5762 void concatwav(int wav[], int **dh, int
|
else if (ij < modmincovj) |
else if (ij < modmincovj) |
modmincovj=ij; |
modmincovj=ij; |
if (ij <0 || ij >1 ){ |
if (ij <0 || ij >1 ){ |
printf("Information, IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); |
printf("ERROR, IMaCh doesn't treat covariate with missing values V%d=-1, individual %d will be skipped.\n",Tvar[k],i); |
fprintf(ficlog,"Information, currently IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i); |
fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=-1, individual %d will be skipped.\n",Tvar[k],i); |
|
fflush(ficlog); |
|
exit(1); |
} |
} |
if ((ij < -1) || (ij > NCOVMAX)){ |
if ((ij < -1) || (ij > NCOVMAX)){ |
printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); |
printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX ); |
Line 5552 void concatwav(int wav[], int **dh, int
|
Line 5840 void concatwav(int wav[], int **dh, int
|
break; |
break; |
} /* end switch */ |
} /* end switch */ |
} /* end dummy test */ |
} /* end dummy test */ |
|
if(Dummy[k]==1 && Typevar[k] !=1){ /* Dummy covariate and not age product */ |
|
for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/ |
|
if(isnan(covar[Tvar[k]][i])){ |
|
printf("ERROR, IMaCh doesn't treat fixed quantitative covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i); |
|
fprintf(ficlog,"ERROR, currently IMaCh doesn't treat covariate with missing values V%d=., individual %d will be skipped.\n",Tvar[k],i); |
|
fflush(ficlog); |
|
exit(1); |
|
} |
|
} |
|
} |
} /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/ |
} /* end of loop on model-covariate k. nbcode[Tvark][1]=-1, nbcode[Tvark][1]=0 and nbcode[Tvark][2]=1 sets the value of covariate k*/ |
|
|
for (k=-1; k< maxncov; k++) Ndum[k]=0; |
for (k=-1; k< maxncov; k++) Ndum[k]=0; |
Line 5865 void concatwav(int wav[], int **dh, int
|
Line 6163 void concatwav(int wav[], int **dh, int
|
varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; |
varhe[ij][ji][(int)age] += doldm[ij][ji]*hf*hf; |
} |
} |
} |
} |
|
/* if((int)age ==50){ */ |
|
/* printf(" age=%d cij=%d nres=%d varhe[%d][%d]=%f ",(int)age, cij, nres, 1,2,varhe[1][2]); */ |
|
/* } */ |
/* Computing expectancies */ |
/* Computing expectancies */ |
hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij,nres); |
hpxij(p3matm,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, cij,nres); |
for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
Line 6528 void varprob(char optionfilefiname[], do
|
Line 6828 void varprob(char optionfilefiname[], do
|
int k2, l2, j1, z1; |
int k2, l2, j1, z1; |
int k=0, l; |
int k=0, l; |
int first=1, first1, first2; |
int first=1, first1, first2; |
|
int nres=0; /* New */ |
double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; |
double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; |
double **dnewm,**doldm; |
double **dnewm,**doldm; |
double *xp; |
double *xp; |
Line 6616 To be simple, these graphs help to under
|
Line 6917 To be simple, these graphs help to under
|
if (cptcovn<1) {tj=1;ncodemax[1]=1;} |
if (cptcovn<1) {tj=1;ncodemax[1]=1;} |
j1=0; |
j1=0; |
for(j1=1; j1<=tj;j1++){ /* For each valid combination of covariates or only once*/ |
for(j1=1; j1<=tj;j1++){ /* For each valid combination of covariates or only once*/ |
|
for(nres=1;nres <=1; nres++){ /* For each resultline */ |
|
/* for(nres=1;nres <=nresult; nres++){ /\* For each resultline *\/ */ |
if (cptcovn>0) { |
if (cptcovn>0) { |
fprintf(ficresprob, "\n#********** Variable "); |
fprintf(ficresprob, "\n#********** Variable "); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
Line 6630 To be simple, these graphs help to under
|
Line 6933 To be simple, these graphs help to under
|
|
|
|
|
fprintf(fichtmcov, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable "); |
fprintf(fichtmcov, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable "); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
/* for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); */ |
|
for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtmcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">"); |
fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">"); |
|
|
fprintf(ficresprobcor, "\n#********** Variable "); |
fprintf(ficresprobcor, "\n#********** Variable "); |
Line 6650 To be simple, these graphs help to under
|
Line 6954 To be simple, these graphs help to under
|
cov[2]=age; |
cov[2]=age; |
if(nagesqr==1) |
if(nagesqr==1) |
cov[3]= age*age; |
cov[3]= age*age; |
for (k=1; k<=cptcovn;k++) { |
/* for (k=1; k<=cptcovn;k++) { */ |
cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; |
/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,k)]; */ |
|
for (k=1; k<=nsd;k++) { /* For single dummy covariates only */ |
|
/* Here comes the value of the covariate 'j1' after renumbering k with single dummy covariates */ |
|
cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(j1,k)]; |
/*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4 |
/*cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];*//* j1 1 2 3 4 |
* 1 1 1 1 1 |
* 1 1 1 1 1 |
* 2 2 1 1 1 |
* 2 2 1 1 1 |
Line 6659 To be simple, these graphs help to under
|
Line 6966 To be simple, these graphs help to under
|
*/ |
*/ |
/* nbcode[1][1]=0 nbcode[1][2]=1;*/ |
/* nbcode[1][1]=0 nbcode[1][2]=1;*/ |
} |
} |
/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
/* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */ |
for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; |
/* ) p nbcode[Tvar[Tage[k]]][(1 & (ij-1) >> (k-1))+1] */ |
for (k=1; k<=cptcovprod;k++) |
/*for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; |
for (k=1; k<=cptcovage;k++){ /* For product with age */ |
|
if(Dummy[Tage[k]]==2){ /* dummy with age */ |
|
cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(j1,k)]*cov[2]; |
|
/* cov[++k1]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ |
|
} else if(Dummy[Tage[k]]==3){ /* quantitative with age */ |
|
printf("Internal IMaCh error, don't know which value for quantitative covariate with age, Tage[k]%d, k=%d, Tvar[Tage[k]]=V%d, age=%d\n",Tage[k],k ,Tvar[Tage[k]], (int)cov[2]); |
|
exit(1); |
|
/* cov[2+nagesqr+Tage[k]]=meanq[k]/idq[k]*cov[2];/\* Using the mean of quantitative variable Tvar[Tage[k]] /\* Tqresult[nres][k]; *\/ */ |
|
/* cov[++k1]=Tqresult[nres][k]; */ |
|
} |
|
/* cov[2+Tage[k]+nagesqr]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2]; */ |
|
} |
|
for (k=1; k<=cptcovprod;k++){/* For product without age */ |
|
if(Dummy[Tvard[k][1]==0]){ |
|
if(Dummy[Tvard[k][2]==0]){ |
|
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,k)] * nbcode[Tvard[k][2]][codtabm(j1,k)]; |
|
/* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */ |
|
}else{ /* Should we use the mean of the quantitative variables? */ |
|
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(j1,k)] * Tqresult[nres][k]; |
|
/* cov[++k1]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k]; */ |
|
} |
|
}else{ |
|
if(Dummy[Tvard[k][2]==0]){ |
|
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(j1,k)] * Tqinvresult[nres][Tvard[k][1]]; |
|
/* cov[++k1]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]]; */ |
|
}else{ |
|
cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; |
|
/* cov[++k1]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]]; */ |
|
} |
|
} |
|
/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)]; */ |
|
} |
|
/* For each age and combination of dummy covariates we slightly move the parameters of delti in order to get the gradient*/ |
for(theta=1; theta <=npar; theta++){ |
for(theta=1; theta <=npar; theta++){ |
for(i=1; i<=npar; i++) |
for(i=1; i<=npar; i++) |
xp[i] = x[i] + (i==theta ?delti[theta]:(double)0); |
xp[i] = x[i] + (i==theta ?delti[theta]:(double)0); |
Line 6849 To be simple, these graphs help to under
|
Line 7186 To be simple, these graphs help to under
|
} /* k12 */ |
} /* k12 */ |
} /*l1 */ |
} /*l1 */ |
}/* k1 */ |
}/* k1 */ |
|
} /* loop on nres */ |
} /* loop on combination of covariates j1 */ |
} /* loop on combination of covariates j1 */ |
free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); |
free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); |
free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); |
free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); |
Line 6871 void printinghtml(char fileresu[], char
|
Line 7209 void printinghtml(char fileresu[], char
|
double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \ |
double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \ |
double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){ |
double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){ |
int jj1, k1, i1, cpt, k4, nres; |
int jj1, k1, i1, cpt, k4, nres; |
|
/* In fact some results are already printed in fichtm which is open */ |
fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \ |
fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \ |
<li><a href='#secondorder'>Result files (second order (variance)</a>\n \ |
<li><a href='#secondorder'>Result files (second order (variance)</a>\n \ |
</ul>"); |
</ul>"); |
fprintf(fichtm,"<ul><li> model=1+age+%s\n \ |
/* fprintf(fichtm,"<ul><li> model=1+age+%s\n \ */ |
</ul>", model); |
/* </ul>", model); */ |
fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n"); |
fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n"); |
fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n", |
fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n", |
jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm")); |
jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm")); |
Line 6909 void printinghtml(char fileresu[], char
|
Line 7247 void printinghtml(char fileresu[], char
|
m=pow(2,cptcoveff); |
m=pow(2,cptcoveff); |
if (cptcovn < 1) {m=1;ncodemax[1]=1;} |
if (cptcovn < 1) {m=1;ncodemax[1]=1;} |
|
|
fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>"); |
fprintf(fichtm," \n<ul><li><b>Graphs (first order)</b></li><p>"); |
|
|
jj1=0; |
jj1=0; |
|
|
Line 6944 void printinghtml(char fileresu[], char
|
Line 7282 void printinghtml(char fileresu[], char
|
fprintf(fichtm,"</a></li>"); |
fprintf(fichtm,"</a></li>"); |
} /* cptcovn >0 */ |
} /* cptcovn >0 */ |
} |
} |
fprintf(fichtm," \n</ul>"); |
fprintf(fichtm," \n</ul>"); |
|
|
jj1=0; |
jj1=0; |
|
|
Line 6978 void printinghtml(char fileresu[], char
|
Line 7316 void printinghtml(char fileresu[], char
|
} |
} |
|
|
/* if(nqfveff+nqtveff 0) */ /* Test to be done */ |
/* if(nqfveff+nqtveff 0) */ /* Test to be done */ |
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
fprintf(fichtm," (model=%s) ************\n<hr size=\"2\" color=\"#EC5E5E\">",model); |
if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
fprintf(fichtm,"\n<h3>Combination (%d) ignored because no cases </h3>\n",k1); |
fprintf(fichtm,"\n<h3>Combination (%d) ignored because no cases </h3>\n",k1); |
printf("\nCombination (%d) ignored because no cases \n",k1); |
printf("\nCombination (%d) ignored because no cases \n",k1); |
Line 7023 divided by h: <sub>h</sub>P<sub>ij</sub>
|
Line 7361 divided by h: <sub>h</sub>P<sub>ij</sub>
|
if(prevfcast==1){ |
if(prevfcast==1){ |
/* Projection of prevalence up to period (forward stable) prevalence in each health state */ |
/* Projection of prevalence up to period (forward stable) prevalence in each health state */ |
for(cpt=1; cpt<=nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) forward prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) forward prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a>", dateprev1, dateprev2, mobilavproj, dateprojd, dateprojf, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres); |
<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, dateprojd, dateprojf, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres); |
fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"F_"),subdirf2(optionfilefiname,"F_")); |
|
fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">", |
|
subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres); |
} |
} |
} |
} |
if(prevbcast==1){ |
if(prevbcast==1){ |
Line 7033 divided by h: <sub>h</sub>P<sub>ij</sub>
|
Line 7373 divided by h: <sub>h</sub>P<sub>ij</sub>
|
fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \ |
fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \ |
from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \ |
from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \ |
account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \ |
account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \ |
with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \ |
with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a>", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres); |
<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres); |
fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"FB_"),subdirf2(optionfilefiname,"FB_")); |
|
fprintf(fichtm," <img src=\"%s_%d-%d-%d.svg\">", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres); |
} |
} |
} |
} |
|
|
for(cpt=1; cpt<=nlstate;cpt++) { |
for(cpt=1; cpt<=nlstate;cpt++) { |
fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a> <br> \ |
fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a>",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres); |
<img src=\"%s_%d-%d-%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres); |
fprintf(fichtm," (data from text file <a href=\"%s.txt\"> %s.txt</a>)\n<br>",subdirf2(optionfilefiname,"E_"),subdirf2(optionfilefiname,"E_")); |
|
fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">", subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres ); |
} |
} |
/* } /\* end i1 *\/ */ |
/* } /\* end i1 *\/ */ |
}/* End k1 */ |
}/* End k1 */ |
Line 7092 See page 'Matrix of variance-covariance
|
Line 7434 See page 'Matrix of variance-covariance
|
/* else */ |
/* else */ |
/* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */ |
/* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */ |
fflush(fichtm); |
fflush(fichtm); |
fprintf(fichtm," <ul><li><b>Graphs</b></li><p>"); |
|
|
|
m=pow(2,cptcoveff); |
m=pow(2,cptcoveff); |
if (cptcovn < 1) {m=1;ncodemax[1]=1;} |
if (cptcovn < 1) {m=1;ncodemax[1]=1;} |
|
|
|
fprintf(fichtm," <ul><li><b>Graphs (second order)</b></li><p>"); |
|
|
|
jj1=0; |
|
|
|
fprintf(fichtm," \n<ul>"); |
|
for(nres=1; nres <= nresult; nres++) /* For each resultline */ |
|
for(k1=1; k1<=m;k1++){ /* For each combination of covariate */ |
|
if(m != 1 && TKresult[nres]!= k1) |
|
continue; |
|
jj1++; |
|
if (cptcovn > 0) { |
|
fprintf(fichtm,"\n<li><a size=\"1\" color=\"#EC5E5E\" href=\"#rescovsecond"); |
|
for (cpt=1; cpt<=cptcoveff;cpt++){ |
|
fprintf(fichtm,"_V%d=%d_",Tvresult[nres][cpt],(int)Tresult[nres][cpt]); |
|
} |
|
for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ |
|
fprintf(fichtm,"_V%d=%f_",Tvqresult[nres][k4],Tqresult[nres][k4]); |
|
} |
|
fprintf(fichtm,"\">"); |
|
|
|
/* if(nqfveff+nqtveff 0) */ /* Test to be done */ |
|
fprintf(fichtm,"************ Results for covariates"); |
|
for (cpt=1; cpt<=cptcoveff;cpt++){ |
|
fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],(int)Tresult[nres][cpt]); |
|
} |
|
for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ |
|
fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); |
|
} |
|
if(invalidvarcomb[k1]){ |
|
fprintf(fichtm," Warning Combination (%d) ignored because no cases ",k1); |
|
continue; |
|
} |
|
fprintf(fichtm,"</a></li>"); |
|
} /* cptcovn >0 */ |
|
} |
|
fprintf(fichtm," \n</ul>"); |
|
|
jj1=0; |
jj1=0; |
|
|
for(nres=1; nres <= nresult; nres++){ /* For each resultline */ |
for(nres=1; nres <= nresult; nres++){ /* For each resultline */ |
Line 7106 See page 'Matrix of variance-covariance
|
Line 7484 See page 'Matrix of variance-covariance
|
/* for(i1=1; i1<=ncodemax[k1];i1++){ */ |
/* for(i1=1; i1<=ncodemax[k1];i1++){ */ |
jj1++; |
jj1++; |
if (cptcovn > 0) { |
if (cptcovn > 0) { |
|
fprintf(fichtm,"\n<p><a name=\"rescovsecond"); |
|
for (cpt=1; cpt<=cptcoveff;cpt++){ |
|
fprintf(fichtm,"_V%d=%d_",Tvresult[nres][cpt],(int)Tresult[nres][cpt]); |
|
} |
|
for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ |
|
fprintf(fichtm,"_V%d=%f_",Tvqresult[nres][k4],Tqresult[nres][k4]); |
|
} |
|
fprintf(fichtm,"\"</a>"); |
|
|
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates"); |
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates"); |
for (cpt=1; cpt<=cptcoveff;cpt++) /**< cptcoveff number of variables */ |
for (cpt=1; cpt<=cptcoveff;cpt++){ /**< cptcoveff number of variables */ |
fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]); |
fprintf(fichtm," V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]); |
|
printf(" V%d=%d ",Tvresult[nres][cpt],Tresult[nres][cpt]);fflush(stdout); |
/* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */ |
/* fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); */ |
|
} |
for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ |
for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ |
fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); |
fprintf(fichtm," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); |
} |
} |
|
|
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
fprintf(fichtm," (model=%s) ************\n<hr size=\"2\" color=\"#EC5E5E\">",model); |
|
|
if(invalidvarcomb[k1]){ |
if(invalidvarcomb[k1]){ |
fprintf(fichtm,"\n<h4>Combination (%d) ignored because no cases </h4>\n",k1); |
fprintf(fichtm,"\n<h4>Combination (%d) ignored because no cases </h4>\n",k1); |
Line 7123 See page 'Matrix of variance-covariance
|
Line 7512 See page 'Matrix of variance-covariance
|
} |
} |
for(cpt=1; cpt<=nlstate;cpt++) { |
for(cpt=1; cpt<=nlstate;cpt++) { |
fprintf(fichtm,"\n<br>- Observed (cross-sectional with mov_average=%d) and period (incidence based) \ |
fprintf(fichtm,"\n<br>- Observed (cross-sectional with mov_average=%d) and period (incidence based) \ |
prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d-%d.svg\"> %s_%d-%d-%d.svg</a>\n <br>\ |
prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d-%d.svg\"> %s_%d-%d-%d.svg</a>",mobilav,cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres); |
<img src=\"%s_%d-%d-%d.svg\">",mobilav,cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres); |
fprintf(fichtm," (data from text file <a href=\"%s\">%s</a>)\n <br>",subdirf2(fileresu,"VPL_"),subdirf2(fileresu,"VPL_")); |
|
fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"V_"), cpt,k1,nres); |
} |
} |
fprintf(fichtm,"\n<br>- Total life expectancy by age and \ |
fprintf(fichtm,"\n<br>- Total life expectancy by age and \ |
health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \ |
health expectancies in each live states (1 to %d). If popbased=1 the smooth (due to the model) \ |
true period expectancies (those weighted with period prevalences are also\ |
true period expectancies (those weighted with period prevalences are also\ |
drawn in addition to the population based expectancies computed using\ |
drawn in addition to the population based expectancies computed using\ |
observed and cahotic prevalences: <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>\n<br>\ |
observed and cahotic prevalences: <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>",nlstate, subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres); |
<img src=\"%s_%d-%d.svg\">",subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres); |
fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>) \n<br>",subdirf2(optionfilefiname,"T_"),subdirf2(optionfilefiname,"T_")); |
|
fprintf(fichtm,"<img src=\"%s_%d-%d.svg\">",subdirf2(optionfilefiname,"E_"),k1,nres); |
/* } /\* end i1 *\/ */ |
/* } /\* end i1 *\/ */ |
}/* End k1 */ |
}/* End k1 */ |
}/* End nres */ |
}/* End nres */ |
Line 7248 void printinggnuplot(char fileresu[], ch
|
Line 7639 void printinggnuplot(char fileresu[], ch
|
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres); |
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres); |
fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres); |
fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres); |
/* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */ |
/* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */ |
fprintf(ficgp,"set title \"Alive state %d %s\" font \"Helvetica,12\"\n",cpt,gplotlabel); |
fprintf(ficgp,"set title \"Alive state %d %s model=%s\" font \"Helvetica,12\"\n",cpt,gplotlabel,model); |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); |
/* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */ |
/* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */ |
/* k1-1 error should be nres-1*/ |
/* k1-1 error should be nres-1*/ |
Line 7669 set ter svg size 640, 480\nunset log y\n
|
Line 8060 set ter svg size 640, 480\nunset log y\n
|
fprintf(ficgp,", '' "); |
fprintf(ficgp,", '' "); |
/* l=(nlstate+ndeath)*(i-1)+1; */ |
/* l=(nlstate+ndeath)*(i-1)+1; */ |
l=(nlstate+ndeath)*(cpt-1)+1; /* fixed for i; cpt=1 1, cpt=2 1+ nlstate+ndeath, 1+2*(nlstate+ndeath) */ |
l=(nlstate+ndeath)*(cpt-1)+1; /* fixed for i; cpt=1 1, cpt=2 1+ nlstate+ndeath, 1+2*(nlstate+ndeath) */ |
/* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */ |
/* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */ |
/* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */ |
/* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */ |
fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+i-1); /* To be verified */ |
fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+i-1); /* To be verified */ |
/* for (j=2; j<= nlstate ; j ++) */ |
/* for (j=2; j<= nlstate ; j ++) */ |
/* fprintf(ficgp,"+$%d",k+l+j-1); */ |
/* fprintf(ficgp,"+$%d",k+l+j-1); */ |
Line 8020 set ter svg size 640, 480\nunset log y\n
|
Line 8411 set ter svg size 640, 480\nunset log y\n
|
for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */ |
for(j=1; j <=cptcovt; j++) { /* For each covariate of the simplified model */ |
/* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */ |
/* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */ |
if(cptcovage >0){ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */ |
if(cptcovage >0){ /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */ |
if(j==Tage[ij]) { /* Product by age To be looked at!!*/ |
if(j==Tage[ij]) { /* Product by age To be looked at!!*//* Bug valgrind */ |
if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */ |
if(ij <=cptcovage) { /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, 2 V5 and V1 */ |
if(DummyV[j]==0){ |
if(DummyV[j]==0){/* Bug valgrind */ |
fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);; |
fprintf(ficgp,"+p%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvar[j]]);; |
}else{ /* quantitative */ |
}else{ /* quantitative */ |
fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */ |
fprintf(ficgp,"+p%d*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvar[j]]); /* Tqinvresult in decoderesult */ |
Line 9212 int readdata(char datafile[], int firsto
|
Line 9603 int readdata(char datafile[], int firsto
|
|
|
DummyV=ivector(1,NCOVMAX); /* 1 to 3 */ |
DummyV=ivector(1,NCOVMAX); /* 1 to 3 */ |
FixedV=ivector(1,NCOVMAX); /* 1 to 3 */ |
FixedV=ivector(1,NCOVMAX); /* 1 to 3 */ |
|
for(v=1;v<NCOVMAX;v++){ |
|
DummyV[v]=0; |
|
FixedV[v]=0; |
|
} |
|
|
for(v=1; v <=ncovcol;v++){ |
for(v=1; v <=ncovcol;v++){ |
DummyV[v]=0; |
DummyV[v]=0; |
Line 9255 int readdata(char datafile[], int firsto
|
Line 9650 int readdata(char datafile[], int firsto
|
else if( line[0] == (char)0xFF && line[1] == (char)0xFE) |
else if( line[0] == (char)0xFF && line[1] == (char)0xFE) |
{ |
{ |
noffset=noffset+2; |
noffset=noffset+2; |
printf("# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout); |
printf("# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout); |
fprintf(ficlog,"# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile); |
fprintf(ficlog,"# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile); |
fflush(ficlog); return 1; |
fflush(ficlog); return 1; |
} |
} |
else if( line[0] == 0 && line[1] == 0) |
else if( line[0] == 0 && line[1] == 0) |
{ |
{ |
if( line[2] == (char)0xFE && line[3] == (char)0xFF){ |
if( line[2] == (char)0xFE && line[3] == (char)0xFF){ |
noffset=noffset+4; |
noffset=noffset+4; |
printf("# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout); |
printf("# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile);fflush(stdout); |
fprintf(ficlog,"# Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile); |
fprintf(ficlog,"# Error Data file '%s' is a huge UTF16BE BOM file, please convert to UTF8 or ascii file (for example with dos2unix) and rerun.\n",datafile); |
fflush(ficlog); return 1; |
fflush(ficlog); return 1; |
} |
} |
} else{ |
} else{ |
Line 9352 int readdata(char datafile[], int firsto
|
Line 9747 int readdata(char datafile[], int firsto
|
} |
} |
if(lval <-1 || lval >1){ |
if(lval <-1 || lval >1){ |
printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ |
printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ |
Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ |
Should be a value of %d(nth) covariate of wave %d (0 should be the value for the reference and 1\n \ |
for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ |
for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ |
For example, for multinomial values like 1, 2 and 3,\n \ |
For example, for multinomial values like 1, 2 and 3,\n \ |
build V1=0 V2=0 for the reference value (1),\n \ |
build V1=0 V2=0 for the reference value (1),\n \ |
V1=1 V2=0 for (2) \n \ |
V1=1 V2=0 for (2) \n \ |
and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ |
and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ |
output of IMaCh is often meaningless.\n \ |
output of IMaCh is often meaningless.\n \ |
Exiting.\n",lval,linei, i,line,j); |
Exiting.\n",lval,linei, i,line,iv,j); |
fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ |
fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ |
Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ |
Should be a value of %d(nth) covariate of wave %d (0 should be the value for the reference and 1\n \ |
for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ |
for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ |
For example, for multinomial values like 1, 2 and 3,\n \ |
For example, for multinomial values like 1, 2 and 3,\n \ |
build V1=0 V2=0 for the reference value (1),\n \ |
build V1=0 V2=0 for the reference value (1),\n \ |
V1=1 V2=0 for (2) \n \ |
V1=1 V2=0 for (2) \n \ |
and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ |
and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ |
output of IMaCh is often meaningless.\n \ |
output of IMaCh is often meaningless.\n \ |
Exiting.\n",lval,linei, i,line,j);fflush(ficlog); |
Exiting.\n",lval,linei, i,line,iv,j);fflush(ficlog); |
return 1; |
return 1; |
} |
} |
cotvar[j][iv][i]=(double)(lval); |
cotvar[j][iv][i]=(double)(lval); |
Line 9476 int readdata(char datafile[], int firsto
|
Line 9871 int readdata(char datafile[], int firsto
|
cutv(stra, strb, line, ' '); |
cutv(stra, strb, line, ' '); |
if(strb[0]=='.') { /* Missing value */ |
if(strb[0]=='.') { /* Missing value */ |
lval=-1; |
lval=-1; |
|
coqvar[iv][i]=NAN; |
|
covar[ncovcol+iv][i]=NAN; /* including qvar in standard covar for performance reasons */ |
}else{ |
}else{ |
errno=0; |
errno=0; |
/* what_kind_of_number(strb); */ |
/* what_kind_of_number(strb); */ |
Line 9579 int decoderesult ( char resultline[], in
|
Line 9976 int decoderesult ( char resultline[], in
|
char stra[80], strb[80], strc[80], strd[80],stre[80]; |
char stra[80], strb[80], strc[80], strd[80],stre[80]; |
|
|
removefirstspace(&resultline); |
removefirstspace(&resultline); |
printf("decoderesult:%s\n",resultline); |
|
|
|
if (strstr(resultline,"v") !=0){ |
if (strstr(resultline,"v") !=0){ |
printf("Error. 'v' must be in upper case 'V' result: %s ",resultline); |
printf("Error. 'v' must be in upper case 'V' result: %s ",resultline); |
Line 9594 int decoderesult ( char resultline[], in
|
Line 9990 int decoderesult ( char resultline[], in
|
TKresult[nres]=0; /* Combination for the nresult and the model */ |
TKresult[nres]=0; /* Combination for the nresult and the model */ |
return (0); |
return (0); |
} |
} |
|
|
if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */ |
if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */ |
printf("ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs); |
printf("ERROR: the number of variables in this result line, %d, differs from the number of variables used in the model line, %d.\n",j, cptcovs); |
fprintf(ficlog,"ERROR: the number of variable in the resultline, %d, differs from the number of variable used in the model line, %d.\n",j, cptcovs); |
fprintf(ficlog,"ERROR: the number of variables in the resultline, %d, differs from the number of variables used in the model line, %d.\n",j, cptcovs); |
} |
} |
for(k=1; k<=j;k++){ /* Loop on any covariate of the result line */ |
for(k=1; k<=j;k++){ /* Loop on any covariate of the result line */ |
if(nbocc(resultsav,'=') >1){ |
if(nbocc(resultsav,'=') >1){ |
cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' |
cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' (stra is the rest of the resultline to be analyzed in the next loop *//* resultsav= "V4=1 V5=25.1 V3=0" stra= "V5=25.1 V3=0" strb= "V4=1" */ |
resultsav= V4=1 V5=25.1 V3=0 strb=V3=0 stra= V4=1 V5=25.1 */ |
cutl(strc,strd,strb,'='); /* strb:"V4=1" strc="1" strd="V4" */ |
cutl(strc,strd,strb,'='); /* strb:V4=1 strc=1 strd=V4 */ |
|
}else |
}else |
cutl(strc,strd,resultsav,'='); |
cutl(strc,strd,resultsav,'='); |
Tvalsel[k]=atof(strc); /* 1 */ |
Tvalsel[k]=atof(strc); /* 1 */ /* Tvalsel of k is the float value of the kth covariate appearing in this result line */ |
|
|
cutl(strc,stre,strd,'V'); /* strd='V4' strc=4 stre='V' */; |
cutl(strc,stre,strd,'V'); /* strd='V4' strc=4 stre='V' */; |
Tvarsel[k]=atoi(strc); |
Tvarsel[k]=atoi(strc); /* 4 */ /* Tvarsel is the id of the kth covariate in the result line Tvarsel[1] in "V4=1.." is 4.*/ |
/* Typevarsel[k]=1; /\* 1 for age product *\/ */ |
/* Typevarsel[k]=1; /\* 1 for age product *\/ */ |
/* cptcovsel++; */ |
/* cptcovsel++; */ |
if (nbocc(stra,'=') >0) |
if (nbocc(stra,'=') >0) |
strcpy(resultsav,stra); /* and analyzes it */ |
strcpy(resultsav,stra); /* and analyzes it */ |
} |
} |
/* Checking for missing or useless values in comparison of current model needs */ |
/* Checking for missing or useless values in comparison of current model needs */ |
for(k1=1; k1<= cptcovt ;k1++){ /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
for(k1=1; k1<= cptcovt ;k1++){ /* Loop on model. model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
if(Typevar[k1]==0){ /* Single covariate in model */ |
if(Typevar[k1]==0){ /* Single covariate in model *//*0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */ |
match=0; |
match=0; |
for(k2=1; k2 <=j;k2++){/* result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ |
for(k2=1; k2 <=j;k2++){/* Loop on resultline. In result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ |
if(Tvar[k1]==Tvarsel[k2]) {/* Tvar[1]=5 == Tvarsel[2]=5 */ |
if(Tvar[k1]==Tvarsel[k2]) {/* Tvar is coming from the model, Tvarsel from the result. Tvar[1]=5 == Tvarsel[2]=5 */ |
modelresult[k2]=k1;/* modelresult[2]=1 modelresult[1]=2 modelresult[3]=3 modelresult[6]=4 modelresult[9]=5 */ |
modelresult[k2]=k1;/* modelresult[2]=1 modelresult[1]=2 modelresult[3]=3 modelresult[6]=4 modelresult[9]=5 */ |
match=1; |
match=1; /* modelresult of k2 variable of resultline is identical to k1 variable of the model good */ |
break; |
break; |
} |
} |
} |
} |
if(match == 0){ |
if(match == 0){ |
printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); |
printf("Error in result line: V%d is missing in result: %s according to model=%s\n",k1, resultline, model); |
|
fprintf(ficlog,"Error in result line: V%d is missing in result: %s according to model=%s\n",k1, resultline, model); |
|
return 1; |
} |
} |
} |
} |
} |
} |
/* Checking for missing or useless values in comparison of current model needs */ |
/* Checking for missing or useless values in comparison of current model needs */ |
for(k2=1; k2 <=j;k2++){ /* result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ |
for(k2=1; k2 <=j;k2++){ /* Loop on resultline variables: result line V4=1 V5=24.1 V3=1 V2=8 V1=0 */ |
match=0; |
match=0; |
for(k1=1; k1<= cptcovt ;k1++){ /* model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
for(k1=1; k1<= cptcovt ;k1++){ /* loop on model: model line V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ |
if(Typevar[k1]==0){ /* Single */ |
if(Typevar[k1]==0){ /* Single */ |
if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4 */ |
if(Tvar[k1]==Tvarsel[k2]) { /* Tvar[2]=4 == Tvarsel[1]=4 */ |
resultmodel[k1]=k2; /* resultmodel[2]=1 resultmodel[1]=2 resultmodel[3]=3 resultmodel[6]=4 resultmodel[9]=5 */ |
resultmodel[k1]=k2; /* k2th variable of the model corresponds to k1 variable of the model. resultmodel[2]=1 resultmodel[1]=2 resultmodel[3]=3 resultmodel[6]=4 resultmodel[9]=5 */ |
++match; |
++match; |
} |
} |
} |
} |
} |
} |
if(match == 0){ |
if(match == 0){ |
printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); |
printf("Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); |
|
fprintf(ficlog,"Error in result line: %d value missing; result: %s, model=%s\n",k1, resultline, model); |
|
return 1; |
}else if(match > 1){ |
}else if(match > 1){ |
printf("Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model); |
printf("Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model); |
|
fprintf(ficlog,"Error in result line: %d doubled; result: %s, model=%s\n",k2, resultline, model); |
|
return 1; |
} |
} |
} |
} |
|
|
Line 9667 int decoderesult ( char resultline[], in
|
Line 10067 int decoderesult ( char resultline[], in
|
/* V(Tvqresult)=Tqresult V5=25.1 V2=8 Tqresult[nres=1][1]=25.1 */ |
/* V(Tvqresult)=Tqresult V5=25.1 V2=8 Tqresult[nres=1][1]=25.1 */ |
/* V5*age V5 known which value for nres? */ |
/* V5*age V5 known which value for nres? */ |
/* Tqinvresult[2]=8 Tqinvresult[1]=25.1 */ |
/* Tqinvresult[2]=8 Tqinvresult[1]=25.1 */ |
for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* model line */ |
for(k1=1, k=0, k4=0, k4q=0; k1 <=cptcovt;k1++){ /* loop on model line */ |
if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy */ |
if( Dummy[k1]==0 && Typevar[k1]==0 ){ /* Single dummy */ |
k3= resultmodel[k1]; /* resultmodel[2(V4)] = 1=k3 */ |
k3= resultmodel[k1]; /* resultmodel[2(V4)] = 1=k3 */ |
k2=(int)Tvarsel[k3]; /* Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */ |
k2=(int)Tvarsel[k3]; /* Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */ |
Line 9678 int decoderesult ( char resultline[], in
|
Line 10078 int decoderesult ( char resultline[], in
|
printf("Decoderesult Dummy k=%d, V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k2, k3, (int)Tvalsel[k3], k4); |
printf("Decoderesult Dummy k=%d, V(k2=V%d)= Tvalsel[%d]=%d, 2**(%d)\n",k, k2, k3, (int)Tvalsel[k3], k4); |
k4++;; |
k4++;; |
} else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative */ |
} else if( Dummy[k1]==1 && Typevar[k1]==0 ){ /* Single quantitative */ |
k3q= resultmodel[k1]; /* resultmodel[2] = 1=k3 */ |
k3q= resultmodel[k1]; /* resultmodel[1(V5)] = 25.1=k3q */ |
k2q=(int)Tvarsel[k3q]; /* Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 */ |
k2q=(int)Tvarsel[k3q]; /* Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */ |
Tqresult[nres][k4q+1]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */ |
Tqresult[nres][k4q+1]=Tvalsel[k3q]; /* Tqresult[nres][1]=25.1 */ |
Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */ |
Tvqresult[nres][k4q+1]=(int)Tvarsel[k3q]; /* Tvqresult[nres][1]=5 */ |
Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */ |
Tqinvresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* Tqinvresult[nres][5]=25.1 */ |
Line 9702 int decodemodel( char model[], int lasto
|
Line 10102 int decodemodel( char model[], int lasto
|
* - cptcovs number of simple covariates |
* - cptcovs number of simple covariates |
* - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10 |
* - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10 |
* which is a new column after the 9 (ncovcol) variables. |
* which is a new column after the 9 (ncovcol) variables. |
* - if k is a product Vn*Vm covar[k][i] is filled with correct values for each individual |
* - if k is a product Vn*Vm, covar[k][i] is filled with correct values for each individual |
* - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage |
* - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage |
* Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6. |
* Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6. |
* - Tvard[k] p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 . |
* - Tvard[k] p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 . |
*/ |
*/ |
|
/* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */ |
{ |
{ |
int i, j, k, ks, v; |
int i, j, k, ks, v; |
int j1, k1, k2, k3, k4; |
int j1, k1, k2, k3, k4; |
Line 9784 int decodemodel( char model[], int lasto
|
Line 10185 int decodemodel( char model[], int lasto
|
* Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 d1 d1 d2 d2 |
* Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 d1 d1 d2 d2 |
* k= 1 2 3 4 5 6 7 8 9 10 11 12 |
* k= 1 2 3 4 5 6 7 8 9 10 11 12 |
* Tvar[k]= 2 1 3 3 10 11 8 8 5 6 7 8 |
* Tvar[k]= 2 1 3 3 10 11 8 8 5 6 7 8 |
* p Tvar[1]@12={2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} |
* p Tvar[1]@12={2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} |
* p Tprod[1]@2={ 6, 5} |
* p Tprod[1]@2={ 6, 5} |
*p Tvard[1][1]@4= {7, 8, 5, 6} |
*p Tvard[1][1]@4= {7, 8, 5, 6} |
* covar[k][i]= V2 V1 ? V3 V5*V6? V7*V8? ? V8 |
* covar[k][i]= V2 V1 ? V3 V5*V6? V7*V8? ? V8 |
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
*How to reorganize? |
*How to reorganize? Tvars(orted) |
* Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age |
* Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age |
* Tvars {2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} |
* Tvars {2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} |
* {2, 1, 4, 8, 5, 6, 3, 7} |
* {2, 1, 4, 8, 5, 6, 3, 7} |
Line 9814 int decodemodel( char model[], int lasto
|
Line 10215 int decodemodel( char model[], int lasto
|
Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0; |
Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0; |
} |
} |
cptcovage=0; |
cptcovage=0; |
for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */ |
for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model line */ |
cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' |
cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' cutl from left to right |
modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ |
modelsav==V2+V1+V5*age+V4+V3*age strb=V3*age stra=V2+V1V5*age+V4 */ /* <model> "V5+V4+V3+V4*V3+V5*age+V1*age+V1" strb="V5" stra="V4+V3+V4*V3+V5*age+V1*age+V1" */ |
if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */ |
if (nbocc(modelsav,'+')==0) |
|
strcpy(strb,modelsav); /* and analyzes it */ |
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ |
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ |
/*scanf("%d",i);*/ |
/*scanf("%d",i);*/ |
if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */ |
if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V5*age+ V4+V3*age strb=V3*age */ |
cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ |
cutl(strc,strd,strb,'*'); /**< k=1 strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ |
if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */ |
if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */ |
/* covar is not filled and then is empty */ |
/* covar is not filled and then is empty */ |
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+V5*age+V4+V3*age Tvar[5]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */ |
Typevar[k]=1; /* 1 for age product */ |
Typevar[k]=1; /* 1 for age product */ |
cptcovage++; /* Sums the number of covariates which include age as a product */ |
cptcovage++; /* Counts 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; /* V2+V1+V4+V3*age Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */ |
/*printf("stre=%s ", stre);*/ |
/*printf("stre=%s ", stre);*/ |
} else if (strcmp(strd,"age")==0) { /* or age*Vn */ |
} else if (strcmp(strd,"age")==0) { /* or age*Vn */ |
cptcovprod--; |
cptcovprod--; |
Line 9846 int decodemodel( char model[], int lasto
|
Line 10248 int decodemodel( char model[], int lasto
|
Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but |
Tvar[k]=ncovcol+nqv+ntv+nqtv+k1; /* For model-covariate k tells which data-covariate to use but |
because this model-covariate is a construction we invent a new column |
because this model-covariate is a construction we invent a new column |
which is after existing variables ncovcol+nqv+ntv+nqtv + k1 |
which is after existing variables ncovcol+nqv+ntv+nqtv + k1 |
If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2 |
If already ncovcol=4 and model=V2 + V1 +V1*V4 +age*V3 +V3*V2 |
Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */ |
thus after V4 we invent V5 and V6 because age*V3 will be computed in 4 |
|
Tvar[3=V1*V4]=4+1=5 Tvar[5=V3*V2]=4 + 2= 6, Tvar[4=age*V3]=4 etc */ |
Typevar[k]=2; /* 2 for double fixed dummy covariates */ |
Typevar[k]=2; /* 2 for double fixed dummy covariates */ |
cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ |
cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ |
Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ |
Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ |
Tposprod[k]=k1; /* Tpsprod[3]=1, Tposprod[2]=5 */ |
Tposprod[k]=k1; /* Tposprod[3]=1, Tposprod[2]=5 */ |
Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ |
Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ |
Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ |
Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ |
k2=k2+2; /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */ |
k2=k2+2; /* k2 is initialize to -1, We want to store the n and m in Vn*Vm at the end of Tvar */ |
Line 9866 int decodemodel( char model[], int lasto
|
Line 10269 int decodemodel( char model[], int lasto
|
} |
} |
} /* End age is not in the model */ |
} /* End age is not in the model */ |
} /* End if model includes a product */ |
} /* End if model includes a product */ |
else { /* no more sum */ |
else { /* not a product */ |
/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ |
/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ |
/* scanf("%d",i);*/ |
/* scanf("%d",i);*/ |
cutl(strd,strc,strb,'V'); |
cutl(strd,strc,strb,'V'); |
Line 9897 int decodemodel( char model[], int lasto
|
Line 10300 int decodemodel( char model[], int lasto
|
model= V5 + V4 +V3 + V4*V3 + V5*age + V2 + V1*V2 + V1*age + V5*age, 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 7 8 9 |
k = 1 2 3 4 5 6 7 8 9 |
Tvar[k]= 5 4 3 1+1+2+1+1=6 5 2 7 1 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 2 1 1 |
Typevar[k]= 0 0 0 2 1 0 2 1 0 |
Fixed[k] 1 1 1 1 3 0 0 or 2 2 3 |
Fixed[k] 1 1 1 1 3 0 0 or 2 2 3 |
Dummy[k] 1 0 0 0 3 1 1 2 3 |
Dummy[k] 1 0 0 0 3 1 1 2 3 |
Tmodelind[combination of covar]=k; |
Tmodelind[combination of covar]=k; |
Line 9906 int decodemodel( char model[], int lasto
|
Line 10309 int decodemodel( char model[], int lasto
|
/* If Tvar[k] >ncovcol it is a product */ |
/* 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 */ |
/* Computing effective variables, ie used by the model, that is from the cptcovt variables */ |
printf("Model=%s\n\ |
printf("Model=1+age+%s\n\ |
Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ |
Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ |
Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\ |
Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\ |
Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); |
Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); |
fprintf(ficlog,"Model=%s\n\ |
fprintf(ficlog,"Model=1+age+%s\n\ |
Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ |
Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\ |
Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\ |
Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\ |
Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); |
Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); |
Line 9977 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
Line 10380 Dummy[k] 0=dummy (0 1), 1 quantitative (
|
modell[k].subtype= VQ; |
modell[k].subtype= VQ; |
ncovv++; /* Only simple time varying variables */ |
ncovv++; /* Only simple time varying variables */ |
nsq++; |
nsq++; |
TvarsQ[nsq]=Tvar[k]; |
TvarsQ[nsq]=Tvar[k]; /* k=1 Tvar=5 nsq=1 TvarsQ[1]=5 */ |
TvarsQind[nsq]=k; |
TvarsQind[nsq]=k; |
TvarV[ncovv]=Tvar[k]; |
TvarV[ncovv]=Tvar[k]; |
TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */ |
TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Any time varying singele */ |
Line 10843 int hPijx(double *p, int bage, int fage)
|
Line 11246 int hPijx(double *p, int bage, int fage)
|
|
|
/* oldm=oldms;savm=savms; */ |
/* oldm=oldms;savm=savms; */ |
/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */ |
/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); */ |
hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k, nres); |
hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k, nres);/* Bug valgrind */ |
/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */ |
/* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */ |
fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j="); |
fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j="); |
for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
Line 10856 int hPijx(double *p, int bage, int fage)
|
Line 11259 int hPijx(double *p, int bage, int fage)
|
/* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */ |
/* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */ |
for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
for(j=1; j<=nlstate+ndeath;j++) |
for(j=1; j<=nlstate+ndeath;j++) |
fprintf(ficrespijb," %.5f", p3mat[i][j][h]); |
fprintf(ficrespijb," %.5f", p3mat[i][j][h]);/* Bug valgrind */ |
fprintf(ficrespijb,"\n"); |
fprintf(ficrespijb,"\n"); |
} |
} |
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
Line 10909 int main(int argc, char *argv[])
|
Line 11312 int main(int argc, char *argv[])
|
double dum=0.; /* Dummy variable */ |
double dum=0.; /* Dummy variable */ |
double ***p3mat; |
double ***p3mat; |
/* double ***mobaverage; */ |
/* double ***mobaverage; */ |
|
double wald; |
|
|
char line[MAXLINE]; |
char line[MAXLINE]; |
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE]; |
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE]; |
Line 10945 int main(int argc, char *argv[])
|
Line 11349 int main(int argc, char *argv[])
|
double ftolpl=FTOL; |
double ftolpl=FTOL; |
double **prlim; |
double **prlim; |
double **bprlim; |
double **bprlim; |
double ***param; /* Matrix of parameters */ |
double ***param; /* Matrix of parameters, param[i][j][k] param=ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel) |
|
state of origin, state of destination including death, for each covariate: constante, age, and V1 V2 etc. */ |
double ***paramstart; /* Matrix of starting parameter values */ |
double ***paramstart; /* Matrix of starting parameter values */ |
double *p, *pstart; /* p=param[1][1] pstart is for starting values guessed by freqsummary */ |
double *p, *pstart; /* p=param[1][1] pstart is for starting values guessed by freqsummary */ |
double **matcov; /* Matrix of covariance */ |
double **matcov; /* Matrix of covariance */ |
Line 11539 Please run with mle=-1 to get a correct
|
Line 11944 Please run with mle=-1 to get a correct
|
} |
} |
mint=matrix(1,maxwav,firstobs,lastobs); |
mint=matrix(1,maxwav,firstobs,lastobs); |
anint=matrix(1,maxwav,firstobs,lastobs); |
anint=matrix(1,maxwav,firstobs,lastobs); |
s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */ |
s=imatrix(1,maxwav+1,firstobs,lastobs); /* s[i][j] health state for wave i and individual j */ |
|
printf("BUG ncovmodel=%d NCOVMAX=%d 2**ncovmodel=%f BUG\n",ncovmodel,NCOVMAX,pow(2,ncovmodel)); |
tab=ivector(1,NCOVMAX); |
tab=ivector(1,NCOVMAX); |
ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
Line 11603 Please run with mle=-1 to get a correct
|
Line 12009 Please run with mle=-1 to get a correct
|
Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age |
Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age |
4 covariates (3 plus signs) |
4 covariates (3 plus signs) |
Tage[1=V3*age]= 4; Tage[2=age*V4] = 3 |
Tage[1=V3*age]= 4; Tage[2=age*V4] = 3 |
*/ |
*/ |
|
for(i=1;i<NCOVMAX;i++) |
|
Tage[i]=0; |
Tmodelind=ivector(1,NCOVMAX);/** gives the k model position of an |
Tmodelind=ivector(1,NCOVMAX);/** gives the k model position of an |
* individual dummy, fixed or varying: |
* individual dummy, fixed or varying: |
* Tmodelind[Tvaraff[3]]=9,Tvaraff[1]@9={4, |
* Tmodelind[Tvaraff[3]]=9,Tvaraff[1]@9={4, |
Line 11816 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 12224 Title=%s <br>Datafile=%s Firstpass=%d La
|
optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model); |
} |
} |
|
|
fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C) 2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-2013-2016-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br> \ |
fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C) 2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-2013-2016-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br> \ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
<font size=\"2\">IMaCh-%s <br> %s</font> \ |
<font size=\"2\">IMaCh-%s <br> %s</font> \ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
<hr size=\"2\" color=\"#EC5E5E\"> \n\ |
Line 11872 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 12280 Title=%s <br>Datafile=%s Firstpass=%d La
|
<img src=\"%s_.svg\">", subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_")); |
<img src=\"%s_.svg\">", subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_")); |
|
|
|
|
fprintf(fichtm,"\n<h4>Some descriptive statistics </h4>\n<br>Total number of observations=%d <br>\n\ |
fprintf(fichtm,"\n<h4>Some descriptive statistics </h4>\n<br>Number of (used) observations=%d <br>\n\ |
Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n\ |
Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n\ |
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\ |
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\ |
imx,agemin,agemax,jmin,jmax,jmean); |
imx,agemin,agemax,jmin,jmax,jmean); |
Line 12152 Please run with mle=-1 to get a correct
|
Line 12560 Please run with mle=-1 to get a correct
|
|
|
|
|
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); /* Printing model equation */ |
fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
|
|
|
printf("#model= 1 + age "); |
|
fprintf(ficres,"#model= 1 + age "); |
|
fprintf(ficlog,"#model= 1 + age "); |
|
fprintf(fichtm,"\n<ul><li> model=1+age+%s\n \ |
|
</ul>", model); |
|
|
|
fprintf(fichtm,"\n<table style=\"text-align:center; border: 1px solid\">\n"); |
|
fprintf(fichtm, "<tr><th>Model=</th><th>1</th><th>+ age</th>"); |
|
if(nagesqr==1){ |
|
printf(" + age*age "); |
|
fprintf(ficres," + age*age "); |
|
fprintf(ficlog," + age*age "); |
|
fprintf(fichtm, "<th>+ age*age</th>"); |
|
} |
|
for(j=1;j <=ncovmodel-2;j++){ |
|
if(Typevar[j]==0) { |
|
printf(" + V%d ",Tvar[j]); |
|
fprintf(ficres," + V%d ",Tvar[j]); |
|
fprintf(ficlog," + V%d ",Tvar[j]); |
|
fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]); |
|
}else if(Typevar[j]==1) { |
|
printf(" + V%d*age ",Tvar[j]); |
|
fprintf(ficres," + V%d*age ",Tvar[j]); |
|
fprintf(ficlog," + V%d*age ",Tvar[j]); |
|
fprintf(fichtm, "<th>+ V%d*age</th>",Tvar[j]); |
|
}else if(Typevar[j]==2) { |
|
printf(" + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); |
|
fprintf(ficres," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); |
|
fprintf(ficlog," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); |
|
fprintf(fichtm, "<th>+ V%d*V%d</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); |
|
} |
|
} |
|
printf("\n"); |
|
fprintf(ficres,"\n"); |
|
fprintf(ficlog,"\n"); |
|
fprintf(fichtm, "</tr>"); |
|
fprintf(fichtm, "\n"); |
|
|
|
|
for(i=1,jk=1; i <=nlstate; i++){ |
for(i=1,jk=1; i <=nlstate; i++){ |
for(k=1; k <=(nlstate+ndeath); k++){ |
for(k=1; k <=(nlstate+ndeath); k++){ |
if (k != i) { |
if (k != i) { |
|
fprintf(fichtm, "<tr>"); |
printf("%d%d ",i,k); |
printf("%d%d ",i,k); |
fprintf(ficlog,"%d%d ",i,k); |
fprintf(ficlog,"%d%d ",i,k); |
fprintf(ficres,"%1d%1d ",i,k); |
fprintf(ficres,"%1d%1d ",i,k); |
|
fprintf(fichtm, "<td>%1d%1d</td>",i,k); |
for(j=1; j <=ncovmodel; j++){ |
for(j=1; j <=ncovmodel; j++){ |
printf("%12.7f ",p[jk]); |
printf("%12.7f ",p[jk]); |
fprintf(ficlog,"%12.7f ",p[jk]); |
fprintf(ficlog,"%12.7f ",p[jk]); |
fprintf(ficres,"%12.7f ",p[jk]); |
fprintf(ficres,"%12.7f ",p[jk]); |
|
fprintf(fichtm, "<td>%12.7f</td>",p[jk]); |
jk++; |
jk++; |
} |
} |
printf("\n"); |
printf("\n"); |
fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
fprintf(ficres,"\n"); |
fprintf(ficres,"\n"); |
|
fprintf(fichtm, "</tr>\n"); |
} |
} |
} |
} |
} |
} |
|
/* fprintf(fichtm,"</tr>\n"); */ |
|
fprintf(fichtm,"</table>\n"); |
|
fprintf(fichtm, "\n"); |
|
|
if(mle != 0){ |
if(mle != 0){ |
/* Computing hessian and covariance matrix only at a peak of the Likelihood, that is after optimization */ |
/* Computing hessian and covariance matrix only at a peak of the Likelihood, that is after optimization */ |
ftolhess=ftol; /* Usually correct */ |
ftolhess=ftol; /* Usually correct */ |
hesscov(matcov, hess, p, npar, delti, ftolhess, func); |
hesscov(matcov, hess, p, npar, delti, ftolhess, func); |
printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n"); |
printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n"); |
fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n"); |
fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n"); |
|
fprintf(fichtm, "\n<p>The Wald test results are output only if the maximimzation of the Likelihood is performed (mle=1)\n</br>Parameters, Wald tests and Wald-based confidence intervals\n</br> W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n</br> And Wald-based confidence intervals plus and minus 1.96 * W \n </br> It might be better to visualize the covariance matrix. See the page '<a href=\"%s\">Matrix of variance-covariance of one-step probabilities and its graphs</a>'.\n</br>",optionfilehtmcov); |
|
fprintf(fichtm,"\n<table style=\"text-align:center; border: 1px solid\">"); |
|
fprintf(fichtm, "\n<tr><th>Model=</th><th>1</th><th>+ age</th>"); |
|
if(nagesqr==1){ |
|
printf(" + age*age "); |
|
fprintf(ficres," + age*age "); |
|
fprintf(ficlog," + age*age "); |
|
fprintf(fichtm, "<th>+ age*age</th>"); |
|
} |
|
for(j=1;j <=ncovmodel-2;j++){ |
|
if(Typevar[j]==0) { |
|
printf(" + V%d ",Tvar[j]); |
|
fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]); |
|
}else if(Typevar[j]==1) { |
|
printf(" + V%d*age ",Tvar[j]); |
|
fprintf(fichtm, "<th>+ V%d*age</th>",Tvar[j]); |
|
}else if(Typevar[j]==2) { |
|
fprintf(fichtm, "<th>+ V%d*V%d</th>",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]); |
|
} |
|
} |
|
fprintf(fichtm, "</tr>\n"); |
|
|
for(i=1,jk=1; i <=nlstate; i++){ |
for(i=1,jk=1; i <=nlstate; i++){ |
for(k=1; k <=(nlstate+ndeath); k++){ |
for(k=1; k <=(nlstate+ndeath); k++){ |
if (k != i) { |
if (k != i) { |
|
fprintf(fichtm, "<tr valign=top>"); |
printf("%d%d ",i,k); |
printf("%d%d ",i,k); |
fprintf(ficlog,"%d%d ",i,k); |
fprintf(ficlog,"%d%d ",i,k); |
|
fprintf(fichtm, "<td>%1d%1d</td>",i,k); |
for(j=1; j <=ncovmodel; j++){ |
for(j=1; j <=ncovmodel; j++){ |
printf("%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); |
wald=p[jk]/sqrt(matcov[jk][jk]); |
fprintf(ficlog,"%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); |
printf("%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); |
|
fprintf(ficlog,"%12.7f(%12.7f) W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk],sqrt(matcov[jk][jk]), p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); |
|
if(fabs(wald) > 1.96){ |
|
fprintf(fichtm, "<td><b>%12.7f</b></br> (%12.7f)</br>",p[jk],sqrt(matcov[jk][jk])); |
|
}else{ |
|
fprintf(fichtm, "<td>%12.7f (%12.7f)</br>",p[jk],sqrt(matcov[jk][jk])); |
|
} |
|
fprintf(fichtm,"W=%8.3f</br>",wald); |
|
fprintf(fichtm,"[%12.7f;%12.7f]</br></td>", p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); |
jk++; |
jk++; |
} |
} |
printf("\n"); |
printf("\n"); |
fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
|
fprintf(fichtm, "</tr>\n"); |
} |
} |
} |
} |
} |
} |
} /* end of hesscov and Wald tests */ |
} /* end of hesscov and Wald tests */ |
|
fprintf(fichtm,"</table>\n"); |
|
|
/* */ |
/* */ |
fprintf(ficres,"# Scales (for hessian or gradient estimation)\n"); |
fprintf(ficres,"# Scales (for hessian or gradient estimation)\n"); |
Line 12405 Please run with mle=-1 to get a correct
|
Line 12895 Please run with mle=-1 to get a correct
|
} |
} |
|
|
/* Results */ |
/* Results */ |
|
endishere=0; |
nresult=0; |
nresult=0; |
|
parameterline=0; |
do{ |
do{ |
if(!fgets(line, MAXLINE, ficpar)){ |
if(!fgets(line, MAXLINE, ficpar)){ |
endishere=1; |
endishere=1; |
parameterline=14; |
parameterline=15; |
}else if (line[0] == '#') { |
}else if (line[0] == '#') { |
/* If line starts with a # it is a comment */ |
/* If line starts with a # it is a comment */ |
numlinepar++; |
numlinepar++; |
Line 12422 Please run with mle=-1 to get a correct
|
Line 12914 Please run with mle=-1 to get a correct
|
parameterline=11; |
parameterline=11; |
else if(sscanf(line,"prevbackcast=%[^\n]\n",modeltemp)) |
else if(sscanf(line,"prevbackcast=%[^\n]\n",modeltemp)) |
parameterline=12; |
parameterline=12; |
else if(sscanf(line,"result:%[^\n]\n",modeltemp)) |
else if(sscanf(line,"result:%[^\n]\n",modeltemp)){ |
parameterline=13; |
parameterline=13; |
|
} |
else{ |
else{ |
parameterline=14; |
parameterline=14; |
} |
} |
switch (parameterline){ |
switch (parameterline){ /* =0 only if only comments */ |
case 11: |
case 11: |
if((num_filled=sscanf(line,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj)) !=EOF && (num_filled == 8)){ |
if((num_filled=sscanf(line,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj)) !=EOF && (num_filled == 8)){ |
fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); |
fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj); |
Line 12440 Please run with mle=-1 to get a correct
|
Line 12933 Please run with mle=-1 to get a correct
|
prvforecast = 1; |
prvforecast = 1; |
} |
} |
else if((num_filled=sscanf(line,"prevforecast=%d yearsfproj=%lf mobil_average=%d\n",&prevfcast,&yrfproj,&mobilavproj)) !=EOF){/* && (num_filled == 3))*/ |
else if((num_filled=sscanf(line,"prevforecast=%d yearsfproj=%lf mobil_average=%d\n",&prevfcast,&yrfproj,&mobilavproj)) !=EOF){/* && (num_filled == 3))*/ |
printf("prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); |
printf("prevforecast=%d yearsfproj=%.2lf mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); |
fprintf(ficlog,"prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); |
fprintf(ficlog,"prevforecast=%d yearsfproj=%.2lf mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); |
fprintf(ficres,"prevforecast=%d yearsfproj=%lf.2 mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); |
fprintf(ficres,"prevforecast=%d yearsfproj=%.2lf mobil_average=%d\n",prevfcast,yrfproj,mobilavproj); |
prvforecast = 2; |
prvforecast = 2; |
} |
} |
else { |
else { |
Line 12463 Please run with mle=-1 to get a correct
|
Line 12956 Please run with mle=-1 to get a correct
|
prvbackcast = 1; |
prvbackcast = 1; |
} |
} |
else if((num_filled=sscanf(line,"prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",&prevbcast,&yrbproj,&mobilavproj)) ==3){/* && (num_filled == 3))*/ |
else if((num_filled=sscanf(line,"prevbackcast=%d yearsbproj=%lf mobil_average=%d\n",&prevbcast,&yrbproj,&mobilavproj)) ==3){/* && (num_filled == 3))*/ |
printf("prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); |
printf("prevbackcast=%d yearsbproj=%.2lf mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); |
fprintf(ficlog,"prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); |
fprintf(ficlog,"prevbackcast=%d yearsbproj=%.2lf mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); |
fprintf(ficres,"prevbackcast=%d yearsbproj=%lf.2 mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); |
fprintf(ficres,"prevbackcast=%d yearsbproj=%.2lf mobil_average=%d\n",prevbcast,yrbproj,mobilavproj); |
prvbackcast = 2; |
prvbackcast = 2; |
} |
} |
else { |
else { |
Line 12475 Please run with mle=-1 to get a correct
|
Line 12968 Please run with mle=-1 to get a correct
|
} |
} |
break; |
break; |
case 13: |
case 13: |
if((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){ |
num_filled=sscanf(line,"result:%[^\n]\n",resultline); |
if (num_filled == 0){ |
nresult++; /* Sum of resultlines */ |
resultline[0]='\0'; |
printf("Result %d: result:%s\n",nresult, resultline); |
printf("Warning %d: no result line! It should be at minimum 'result: V2=0 V1=1 or result:.\n%s\n", num_filled, line); |
if(nresult > MAXRESULTLINESPONE-1){ |
fprintf(ficlog,"Warning %d: no result line! It should be at minimum 'result: V2=0 V1=1 or result:.\n%s\n", num_filled, line); |
printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINESPONE-1,nresult,rfileres); |
break; |
fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINESPONE-1,nresult,rfileres); |
} else if (num_filled != 1){ |
goto end; |
printf("ERROR %d: result line! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",num_filled, line); |
} |
fprintf(ficlog,"ERROR %d: result line! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",num_filled, line); |
if(!decoderesult(resultline, nresult)){ /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */ |
} |
|
nresult++; /* Sum of resultlines */ |
|
printf("Result %d: result=%s\n",nresult, resultline); |
|
if(nresult > MAXRESULTLINES){ |
|
printf("ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult); |
|
fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\n",MAXRESULTLINES,nresult); |
|
goto end; |
|
} |
|
decoderesult(resultline, nresult); /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */ |
|
fprintf(ficparo,"result: %s\n",resultline); |
fprintf(ficparo,"result: %s\n",resultline); |
fprintf(ficres,"result: %s\n",resultline); |
fprintf(ficres,"result: %s\n",resultline); |
fprintf(ficlog,"result: %s\n",resultline); |
fprintf(ficlog,"result: %s\n",resultline); |
break; |
} else |
case 14: |
|
printf("Error: Unknown command '%s'\n",line); |
|
fprintf(ficlog,"Error: Unknown command '%s'\n",line); |
|
if(ncovmodel >=2 && nresult==0 ){ |
|
printf("ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line); |
|
fprintf(ficlog,"ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line); |
|
} |
|
goto end; |
goto end; |
break; |
break; |
default: |
case 14: |
nresult=1; |
printf("Error: Unknown command '%s'\n",line); |
decoderesult(".",nresult ); /* No covariate */ |
fprintf(ficlog,"Error: Unknown command '%s'\n",line); |
|
if(line[0] == ' ' || line[0] == '\n'){ |
|
printf("It should not be an empty line '%s'\n",line); |
|
fprintf(ficlog,"It should not be an empty line '%s'\n",line); |
|
} |
|
if(ncovmodel >=2 && nresult==0 ){ |
|
printf("ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line); |
|
fprintf(ficlog,"ERROR: no result lines! It should be at minimum 'result: V2=0 V1=1 or result:.' %s\n",line); |
} |
} |
|
/* goto end; */ |
|
break; |
|
case 15: |
|
printf("End of resultlines.\n"); |
|
fprintf(ficlog,"End of resultlines.\n"); |
|
break; |
|
default: /* parameterline =0 */ |
|
nresult=1; |
|
decoderesult(".",nresult ); /* No covariate */ |
} /* End switch parameterline */ |
} /* End switch parameterline */ |
}while(endishere==0); /* End do */ |
}while(endishere==0); /* End do */ |
|
|
Line 12781 Please run with mle=-1 to get a correct
|
Line 13274 Please run with mle=-1 to get a correct
|
for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */ |
for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */ |
if(i1 != 1 && TKresult[nres]!= k) |
if(i1 != 1 && TKresult[nres]!= k) |
continue; |
continue; |
printf("\n#****** Result for:"); |
printf("\n# model %s \n#****** Result for:", model); |
fprintf(ficrest,"\n#****** Result for:"); |
fprintf(ficrest,"\n# model %s \n#****** Result for:", model); |
fprintf(ficlog,"\n#****** Result for:"); |
fprintf(ficlog,"\n# model %s \n#****** Result for:", model); |
for(j=1;j<=cptcoveff;j++){ |
for(j=1;j<=cptcoveff;j++){ |
printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |