/* $Id$
$State$
$Log$
+ Revision 1.186 2015/04/23 12:01:52 brouard
+ Summary: V1*age is working now, version 0.98q1
+
+ Some codes had been disabled in order to simplify and Vn*age was
+ working in the optimization phase, ie, giving correct MLE parameters,
+ but, as usual, outputs were not correct and program core dumped.
+
Revision 1.185 2015/03/11 13:26:42 brouard
Summary: Inclusion of compile and links command line for Intel Compiler
end
*/
+/* #define DEBUG */
+/* #define DEBUGBRENT */
#define POWELL /* Instead of NLOPT */
/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
/* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
-int nvar=0, nforce=0; /* Number of variables, number of forces */
+int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */
/* Number of covariates model=V2+V1+ V3*age+V2*V4 */
int cptcovn=0; /**< cptcovn number of covariates added in the model (excepting constant and age and age*product) */
int cptcovt=0; /**< cptcovt number of covariates added in the model (excepting constant and age) */
double *agedc;
double **covar; /**< covar[j,i], value of jth covariate for individual i,
* covar=matrix(0,NCOVMAX,1,n);
- * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; */
+ * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */
double idx;
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */
int *Ndum; /** Freq of modality (tricode */
return s;
}
+/* char *substrchaine(char *out, char *in, char *chain) */
+/* { */
+/* /\* Substract chain 'chain' from 'in', return and output 'out' *\/ */
+/* char *s, *t; */
+/* t=in;s=out; */
+/* while ((*in != *chain) && (*in != '\0')){ */
+/* *out++ = *in++; */
+/* } */
+
+/* /\* *in matches *chain *\/ */
+/* while ((*in++ == *chain++) && (*in != '\0')){ */
+/* printf("*in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */
+/* } */
+/* in--; chain--; */
+/* while ( (*in != '\0')){ */
+/* printf("Bef *in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */
+/* *out++ = *in++; */
+/* printf("Aft *in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */
+/* } */
+/* *out='\0'; */
+/* out=s; */
+/* return out; */
+/* } */
+char *substrchaine(char *out, char *in, char *chain)
+{
+ /* Substract chain 'chain' from 'in', return and output 'out' */
+ /* in="V1+V1*age+age*age+V2", chain="age*age" */
+
+ char *strloc;
+
+ strcpy (out, in);
+ strloc = strstr(out, chain); /* strloc points to out at age*age+V2 */
+ printf("Bef strloc=%s chain=%s out=%s \n", strloc, chain, out);
+ if(strloc != NULL){
+ /* will affect out */ /* strloc+strlenc(chain)=+V2 */ /* Will also work in Unicode */
+ memmove(strloc,strloc+strlen(chain), strlen(strloc+strlen(chain))+1);
+ /* strcpy (strloc, strloc +strlen(chain));*/
+ }
+ printf("Aft strloc=%s chain=%s in=%s out=%s \n", strloc, chain, in, out);
+ return out;
+}
+
+
char *cutl(char *blocc, char *alocc, char *in, char occ)
{
- /* 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')
- gives blocc="abcdef2ghi" and alocc="j".
+ gives blocc="abcdef" and alocc="ghi2j".
If occ is not found blocc is null and alocc is equal to in. Returns blocc
*/
char *s, *t;
}
char *cutv(char *blocc, char *alocc, char *in, char occ)
{
- /* cuts string in into blocc and alocc where blocc ends before last occurence of char 'occ'
+ /* cuts string in into blocc and alocc where blocc ends before LAST occurence of char 'occ'
and alocc starts after last occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2')
gives blocc="abcdef2ghi" and alocc="j".
If occ is not found blocc is null and alocc is equal to in. Returns alocc
/*****************brent *************************/
double brent(double ax, double bx, double cx, double (*f)(double), double tol, double *xmin)
-{
+{
+ /* Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx is
+ * between ax and cx, and f(bx) is less than both f(ax) and f(cx) ), this routine isolates
+ * the minimum to a fractional precision of about tol using Brent’s method. The abscissa of
+ * the minimum is returned as xmin, and the minimum function value is returned as brent , the
+ * returned function value.
+ */
int iter;
double a,b,d,etemp;
double fu=0,fv,fw,fx;
*/
double ulim,u,r,q, dum;
double fu;
-
- *fa=(*func)(*ax);
- *fb=(*func)(*bx);
+
+ double scale=10.;
+ int iterscale=0;
+
+ *fa=(*func)(*ax); /* xta[j]=pcom[j]+(*ax)*xicom[j]; fa=f(xta[j])*/
+ *fb=(*func)(*bx); /* xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) */
+
+
+ /* while(*fb != *fb){ /\* *ax should be ok, reducing distance to *ax *\/ */
+ /* printf("Warning mnbrak *fb = %lf, *bx=%lf *ax=%lf *fa==%lf iter=%d\n",*fb, *bx, *ax, *fa, iterscale++); */
+ /* *bx = *ax - (*ax - *bx)/scale; */
+ /* *fb=(*func)(*bx); /\* xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) *\/ */
+ /* } */
+
if (*fb > *fa) {
SHFT(dum,*ax,*bx,dum)
SHFT(dum,*fb,*fa,dum)
int j;
double xx,xmin,bx,ax;
double fx,fb,fa;
+
+ double scale=10., axs, xxs, xxss; /* Scale added for infinity */
ncom=n;
pcom=vector(1,n);
pcom[j]=p[j];
xicom[j]=xi[j];
}
- ax=0.0;
- xx=1.0;
- mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Find a bracket a,x,b in direction n=xi ie xicom */
- *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Find a minimum P+lambda n in that direction (lambdamin), with TOL between abscisses */
+
+ axs=0.0;
+ xxss=1; /* 1 and using scale */
+ xxs=1;
+ do{
+ ax=0.;
+ xx= xxs;
+ mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Outputs: xtx[j]=pcom[j]+(*xx)*xicom[j]; fx=f(xtx[j]) */
+ /* brackets with inputs ax=0 and xx=1, but points, pcom=p, and directions values, xicom=xi, are sent via f1dim(x) */
+ /* xt[x,j]=pcom[j]+x*xicom[j] f(ax) = f(xt(a,j=1,n)) = f(p(j) + 0 * xi(j)) and f(xx) = f(xt(x, j=1,n)) = f(p(j) + 1 * xi(j)) */
+ /* Outputs: fa=f(p(j)) and fx=f(p(j) + xxs * xi(j) ) and f(bx)= f(p(j)+ bx* xi(j)) */
+ /* Given input ax=axs and xx=xxs, xx might be too far from ax to get a finite f(xx) */
+ /* Searches on line, outputs (ax, xx, bx) such that fx < min(fa and fb) */
+ /* Find a bracket a,x,b in direction n=xi ie xicom, order may change. Scale is [0:xxs*xi[j]] et non plus [0:xi[j]]*/
+ if (fx != fx){
+ xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */
+ printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx);
+ }
+ }while(fx != fx);
+
+ *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/
+ /* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */
+ /* fmin = f(p[j] + xmin * xi[j]) */
+ /* P+lambda n in that direction (lambdamin), with TOL between abscisses */
+ /* f1dim(xmin): for (j=1;j<=ncom;j++) xt[j]=pcom[j]+xmin*xicom[j]; */
#ifdef DEBUG
printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
#endif
+ printf("linmin end ");
for (j=1;j<=n;j++) {
- xi[j] *= xmin;
- p[j] += xi[j];
+ printf(" before xi[%d]=%12.8f", j,xi[j]);
+ xi[j] *= xmin; /* xi rescaled by xmin: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */
+ if(xxs <1.0)
+ printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs );
+ p[j] += xi[j]; /* Parameters values are updated accordingly */
}
+ printf("\n");
free_vector(xicom,1,n);
free_vector(pcom,1,n);
}
for (j=1;j<=n;j++) pt[j]=p[j];
rcurr_time = time(NULL);
for (*iter=1;;++(*iter)) {
- fp=(*fret);
+ fp=(*fret); /* From former iteration or initial value */
ibig=0;
del=0.0;
rlast_time=rcurr_time;
fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
}
}
- for (i=1;i<=n;i++) {
- for (j=1;j<=n;j++) xit[j]=xi[j][i];
+ for (i=1;i<=n;i++) { /* For each direction i */
+ for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */
fptt=(*fret);
#ifdef DEBUG
printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
#endif
- printf("%d",i);fflush(stdout);
+ printf("%d",i);fflush(stdout); /* print direction (parameter) i */
fprintf(ficlog,"%d",i);fflush(ficlog);
- linmin(p,xit,n,fret,func); /* xit[n] has been loaded for direction i */
+ linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input. Outputs are fret(new point p) p is updated and xit rescaled */
if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions
because that direction will be replaced unless the gain del is small
in comparison with the 'probable' gain, mu^2, with the last average direction.
printf("\n");
fprintf(ficlog,"\n");
#endif
- } /* end i */
+ } /* end loop on each direction i */
+ /* 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 and do not produce *fret any more! */
+ /* New value of last point Pn is not computed, P(n-1) */
if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */
#ifdef DEBUG
int k[2],l;
}
if (directest < 0.0) { /* Then we use it for new direction */
#endif
- linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction.*/
+ linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
for (j=1;j<=n;j++) {
xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */
newm=savm;
/* Covariates have to be included here again */
cov[2]=agefin;
-
+ if(nagesqr==1)
+ cov[3]= agefin*agefin;;
for (k=1; k<=cptcovn;k++) {
- cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
/*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/
}
/*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
- for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2];
+ for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2];
for (k=1; k<=cptcovprod;k++) /* Useless */
- cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
/*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
/*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
int i, j, d, h, k;
double **out, cov[NCOVMAX+1];
double **newm;
+ double agexact;
/* Hstepm could be zero and should return the unit matrix */
for (i=1;i<=nlstate+ndeath;i++)
newm=savm;
/* Covariates have to be included here again */
cov[1]=1.;
- cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
+ agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
for (k=1; k<=cptcovn;k++)
- cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
/* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
- cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2];
+ cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2];
for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
- cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
int s1, s2;
double bbh, survp;
long ipmx;
+ double agexact;
/*extern weight */
/* We are differentiating ll according to initial status */
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
to be observed in j being in i according to the model.
*/
for (k=1; k<=cptcovn;k++){ /* Simple and product covariates without age* products */
- cov[2+k]=covar[Tvar[k]][i];
+ cov[2+nagesqr+k]=covar[Tvar[k]][i];
}
/* 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]
}
for(d=0; d<dh[mi][i]; d++){
newm=savm;
- cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; /* Tage[kk] gives the data-covariate associated with age */
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
}
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
} /* end of individual */
} else if(mle==2){
for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
for(mi=1; mi<= wav[i]-1; mi++){
for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){
}
for(d=0; d<=dh[mi][i]; d++){
newm=savm;
- cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
}
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
} /* end of individual */
} else if(mle==3){ /* exponential inter-extrapolation */
for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
for(mi=1; mi<= wav[i]-1; mi++){
for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){
}
for(d=0; d<dh[mi][i]; d++){
newm=savm;
- cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
}
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
} /* end of individual */
}else if (mle==4){ /* ml=4 no inter-extrapolation */
for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
for(mi=1; mi<= wav[i]-1; mi++){
for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){
}
for(d=0; d<dh[mi][i]; d++){
newm=savm;
- cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
}
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
} /* end of individual */
}else{ /* ml=5 no inter-extrapolation no jackson =0.8a */
for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
for(mi=1; mi<= wav[i]-1; mi++){
for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){
}
for(d=0; d<dh[mi][i]; d++){
newm=savm;
- cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
}
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
double llt;
int s1, s2;
double bbh, survp;
+ double agexact;
/*extern weight */
/* We are differentiating ll according to initial status */
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
for(k=1; k<=nlstate; k++) ll[k]=0.;
for (i=1,ipmx=0, sw=0.; i<=imx; i++){
- for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
+ for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
for(mi=1; mi<= wav[i]-1; mi++){
for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){
}
for(d=0; d<dh[mi][i]; d++){
newm=savm;
- cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
+ cov[2]=agexact;
+ if(nagesqr==1)
+ cov[3]= agexact*agexact;
for (kk=1; kk<=cptcovage;kk++) {
- cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
}
+
/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
/* getting the maximum value of the modality of the covariate
(should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
female is 1, then modmaxcovj=1.*/
- }
+ } /* end for loop on individuals */
printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj);
cptcode=modmaxcovj;
/* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */
/*for (i=0; i<=cptcode; i++) {*/
for (i=modmincovj; i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*//* For each value of the modality of model-cov j */
- printf("Frequencies of covariates %d V%d %d\n", j, Tvar[j], Ndum[i]);
+ printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], i, Ndum[i]);
if( Ndum[i] != 0 ){ /* Counts if nobody answered, empty modality */
ncodemax[j]++; /* ncodemax[j]= Number of non-null modalities of the j th covariate. */
}
for (k=-1; k< maxncov; k++) Ndum[k]=0;
- for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */
+ for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */
/* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
ij=Tvar[i]; /* Tvar might be -1 if status was unknown */
- Ndum[ij]++;
+ Ndum[ij]++; /* Might be supersed V1 + V1*age */
}
ij=1;
gm=vector(1,(nlstate)*(nlstate+ndeath));
for (age=bage; age<=fage; age ++){
cov[2]=age;
+ if(nagesqr==1)
+ cov[3]= age*age;
for (k=1; k<=cptcovn;k++) {
- cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4
+ cov[2+nagesqr+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4
* 1 1 1 1 1
* 2 2 1 1 1
* 3 1 2 1 1
/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2];
for (k=1; k<=cptcovprod;k++)
- cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+ cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
for(theta=1; theta <=npar; theta++){
} /* end covariate */
/* proba elementaires */
+ fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n");
for(i=1,jk=1; i <=nlstate; i++){
+ fprintf(ficgp,"# initial state %d\n",i);
for(k=1; k <=(nlstate+ndeath); k++){
if (k != i) {
+ fprintf(ficgp,"# current state %d\n",k);
for(j=1; j <=ncovmodel; j++){
- fprintf(ficgp,"p%d=%f ",jk,p[jk]);
+ fprintf(ficgp,"p%d=%f; ",jk,p[jk]);
jk++;
- fprintf(ficgp,"\n");
}
+ fprintf(ficgp,"\n");
}
}
}
+ fprintf(ficgp,"##############\n#\n");
+
/*goto avoid;*/
+ fprintf(ficgp,"\n##############\n#Graphics of of probabilities or incidences\n#############\n");
+ fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
+ fprintf(ficgp,"# logi(p12/p11)=p1 +p2*age +p3*age*age+ p4*V1+ p5*V1*age\n");
+ fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n");
+ fprintf(ficgp,"# logi(p13/p11)=p6 +p7*age +p8*age*age+ p9*V1+ p10*V1*age\n");
+ fprintf(ficgp,"# p12+p13+p14+p11=1=p11(1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n");
+ fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age)+...)\n");
+ fprintf(ficgp,"# p11=1/(1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n");
+ fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age)+...)\n");
+ fprintf(ficgp,"# p12=exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)/\n");
+ fprintf(ficgp,"# (1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n");
+ fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age))\n");
+ fprintf(ficgp,"# +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n");
+ fprintf(ficgp,"#\n");
for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/
+ fprintf(ficgp,"# ng=%d\n",ng);
+ fprintf(ficgp,"# jk=1 to 2^%d=%d\n",cptcoveff,m);
for(jk=1; jk <=m; jk++) {
+ fprintf(ficgp,"# jk=%d\n",jk);
fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"pe"),jk,ng);
if (ng==2)
fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n");
for(k=1; k<=(nlstate+ndeath); k++) {
if (k != k2){
if(ng==2)
- fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
+ if(nagesqr==0)
+ fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1);
+ else /* nagesqr =1 */
+ fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr);
else
- fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
+ if(nagesqr==0)
+ fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
+ else /* nagesqr =1 */
+ fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr);
ij=1;/* To be checked else nbcode[0][0] wrong */
- for(j=3; j <=ncovmodel; j++) {
+ for(j=3; j <=ncovmodel-nagesqr; j++) {
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */
- fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+ fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
ij++;
}
else
- fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+ fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
}
fprintf(ficgp,")/(1");
- for(k1=1; k1 <=nlstate; k1++){
- fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+ for(k1=1; k1 <=nlstate; k1++){
+ if(nagesqr==0)
+ fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
+ else /* nagesqr =1 */
+ fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
+
ij=1;
- for(j=3; j <=ncovmodel; j++){
+ for(j=3; j <=ncovmodel-nagesqr; j++){
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
- fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+ fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
ij++;
}
else
- fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+ fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
}
fprintf(ficgp,")");
}
}
int decodemodel ( char model[], int lastobs) /**< This routine decode the model and returns:
- * Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age
- * - cptcovt total number of covariates of the model nbocc(+)+1 = 8
- * - cptcovn or number of covariates k of the models excluding age*products =6
+ * Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age
+ * - nagesqr = 1 if age*age in the model, otherwise 0.
+ * - cptcovt total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age
+ * - cptcovn or number of covariates k of the models excluding age*products =6 and age*age
* - cptcovage number of covariates with age*products =2
* - 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
int j1, k1, k2;
char modelsav[80];
char stra[80], strb[80], strc[80], strd[80],stre[80];
+ char *strpt;
/*removespace(model);*/
if (strlen(model) >1){ /* If there is at least 1 covariate */
j=0, j1=0, k1=0, k2=-1, ks=0, cptcovn=0;
- j=nbocc(model,'+'); /**< j=Number of '+' */
- j1=nbocc(model,'*'); /**< j1=Number of '*' */
- cptcovs=j+1-j1; /**< Number of simple covariates V1+V2*age+V3 +V3*V4=> V1 + V3 =2 */
- cptcovt= j+1; /* Number of total covariates in the model V1 + V2*age+ V3 + V3*V4=> 4*/
- /* including age products which are counted in cptcovage.
- * but the covariates which are products must be treated separately: ncovn=4- 2=2 (V1+V3). */
- cptcovprod=j1; /**< Number of products V1*V2 +v3*age = 2 */
- cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */
- strcpy(modelsav,model);
if (strstr(model,"AGE") !=0){
- printf("Error. AGE must be in lower case 'age' model=%s ",model);
- fprintf(ficlog,"Error. AGE must be in lower case model=%s ",model);fflush(ficlog);
+ printf("Error. AGE must be in lower case 'age' model=1+age+%s ",model);
+ fprintf(ficlog,"Error. AGE must be in lower case model=1+age+%s ",model);fflush(ficlog);
return 1;
}
if (strstr(model,"v") !=0){
fprintf(ficlog,"Error. 'v' must be in upper case model=%s ",model);fflush(ficlog);
return 1;
}
+ strcpy(modelsav,model);
+ if ((strpt=strstr(model,"age*age")) !=0){
+ printf(" strpt=%s, model=%s\n",strpt, model);
+ if(strpt != model){
+ printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+ 'model=1+age+age*age+V1' or 'model=1+age+age*age+V1+V1*age', please swap as well as \n \
+ corresponding column of parameters.\n",model);
+ fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
+ 'model=1+age+age*age+V1' or 'model=1+age+age*age+V1+V1*age', please swap as well as \n \
+ corresponding column of parameters.\n",model); fflush(ficlog);
+ return 1;
+ }
+
+ nagesqr=1;
+ if (strstr(model,"+age*age") !=0)
+ substrchaine(modelsav, model, "+age*age");
+ else if (strstr(model,"age*age+") !=0)
+ substrchaine(modelsav, model, "age*age+");
+ else
+ substrchaine(modelsav, model, "age*age");
+ }else
+ nagesqr=0;
+ if (strlen(modelsav) >1){
+ j=nbocc(modelsav,'+'); /**< j=Number of '+' */
+ j1=nbocc(modelsav,'*'); /**< j1=Number of '*' */
+ cptcovs=j+1-j1; /**< Number of simple covariates V1+V1*age+V3 +V3*V4+age*age=> V1 + V3 =2 */
+ cptcovt= j+1; /* Number of total covariates in the model, not including
+ * cst, age and age*age
+ * V1+V1*age+ V3 + V3*V4+age*age=> 4*/
+ /* including age products which are counted in cptcovage.
+ * but the covariates which are products must be treated
+ * separately: ncovn=4- 2=2 (V1+V3). */
+ cptcovprod=j1; /**< Number of products V1*V2 +v3*age = 2 */
+ cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */
+
- /* Design
- * V1 V2 V3 V4 V5 V6 V7 V8 V9 Weight
- * < ncovcol=8 >
- * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8
- * k= 1 2 3 4 5 6 7 8
- * cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8
- * covar[k,i], value of kth covariate if not including age for individual i:
- * covar[1][i]= (V2), covar[4][i]=(V3), covar[8][i]=(V8)
- * Tvar[k] # of the kth covariate: Tvar[1]=2 Tvar[4]=3 Tvar[8]=8
- * if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and
- * Tage[++cptcovage]=k
- * if products, new covar are created after ncovcol with k1
- * Tvar[k]=ncovcol+k1; # of the kth covariate product: Tvar[5]=ncovcol+1=10 Tvar[6]=ncovcol+1=11
- * Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product
- * Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8
- * Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2];
- * Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted
- * V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11
- * < ncovcol=8 >
- * 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
- * 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 Tprod[1]@2={ 6, 5}
- *p Tvard[1][1]@4= {7, 8, 5, 6}
- * covar[k][i]= V2 V1 ? V3 V5*V6? V7*V8? ? V8
- * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
- *How to reorganize?
- * 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}
- * {2, 1, 4, 8, 5, 6, 3, 7}
- * Struct []
- */
+ /* Design
+ * V1 V2 V3 V4 V5 V6 V7 V8 V9 Weight
+ * < ncovcol=8 >
+ * Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8
+ * k= 1 2 3 4 5 6 7 8
+ * cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8
+ * covar[k,i], value of kth covariate if not including age for individual i:
+ * covar[1][i]= (V2), covar[4][i]=(V3), covar[8][i]=(V8)
+ * Tvar[k] # of the kth covariate: Tvar[1]=2 Tvar[4]=3 Tvar[8]=8
+ * if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and
+ * Tage[++cptcovage]=k
+ * if products, new covar are created after ncovcol with k1
+ * Tvar[k]=ncovcol+k1; # of the kth covariate product: Tvar[5]=ncovcol+1=10 Tvar[6]=ncovcol+1=11
+ * Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product
+ * Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8
+ * Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2];
+ * Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted
+ * V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11
+ * < ncovcol=8 >
+ * 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
+ * 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 Tprod[1]@2={ 6, 5}
+ *p Tvard[1][1]@4= {7, 8, 5, 6}
+ * covar[k][i]= V2 V1 ? V3 V5*V6? V7*V8? ? V8
+ * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
+ *How to reorganize?
+ * 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}
+ * {2, 1, 4, 8, 5, 6, 3, 7}
+ * Struct []
+ */
- /* This loop fills the array Tvar from the string 'model'.*/
- /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */
- /* modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4 */
- /* k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tage[cptcovage=1]=4 */
- /* k=3 V4 Tvar[k=3]= 4 (from V4) */
- /* k=2 V1 Tvar[k=2]= 1 (from V1) */
- /* k=1 Tvar[1]=2 (from V2) */
- /* k=5 Tvar[5] */
- /* for (k=1; k<=cptcovn;k++) { */
- /* cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */
- /* } */
- /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; */
- /*
- * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */
- for(k=cptcovt; k>=1;k--) /**< Number of covariates */
+ /* This loop fills the array Tvar from the string 'model'.*/
+ /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */
+ /* modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4 */
+ /* k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tage[cptcovage=1]=4 */
+ /* k=3 V4 Tvar[k=3]= 4 (from V4) */
+ /* k=2 V1 Tvar[k=2]= 1 (from V1) */
+ /* k=1 Tvar[1]=2 (from V2) */
+ /* k=5 Tvar[5] */
+ /* for (k=1; k<=cptcovn;k++) { */
+ /* cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */
+ /* } */
+ /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; */
+ /*
+ * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */
+ for(k=cptcovt; k>=1;k--) /**< Number of covariates */
Tvar[k]=0;
- cptcovage=0;
- for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */
- cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'
- modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */
- if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
- /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
- /*scanf("%d",i);*/
- if (strchr(strb,'*')) { /**< Model includes a product V2+V1+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 */
- if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
- /* covar is not filled and then is empty */
- cptcovprod--;
- 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 */
- cptcovage++; /* Sums the number of covariates which include age as a product */
- Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
- /*printf("stre=%s ", stre);*/
- } else if (strcmp(strd,"age")==0) { /* or age*Vn */
- cptcovprod--;
- cutl(stre,strb,strc,'V');
- Tvar[k]=atoi(stre);
- cptcovage++;
- Tage[cptcovage]=k;
- } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/
- /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
+ cptcovage=0;
+ for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */
+ cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'
+ modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */
+ if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
+ /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
+ /*scanf("%d",i);*/
+ if (strchr(strb,'*')) { /**< Model includes a product V2+V1+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 */
+ if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
+ /* covar is not filled and then is empty */
+ cptcovprod--;
+ 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 */
+ cptcovage++; /* Sums the number of covariates which include age as a product */
+ Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
+ /*printf("stre=%s ", stre);*/
+ } else if (strcmp(strd,"age")==0) { /* or age*Vn */
+ cptcovprod--;
+ cutl(stre,strb,strc,'V');
+ Tvar[k]=atoi(stre);
+ cptcovage++;
+ Tage[cptcovage]=k;
+ } else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/
+ /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
+ cptcovn++;
+ cptcovprodnoage++;k1++;
+ cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
+ Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but
+ because this model-covariate is a construction we invent a new column
+ ncovcol + k1
+ 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 */
+ 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 */
+ Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
+ Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
+ k2=k2+2;
+ Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */
+ Tvar[cptcovt+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */
+ for (i=1; i<=lastobs;i++){
+ /* Computes the new covariate which is a product of
+ covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
+ covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
+ }
+ } /* End age is not in the model */
+ } /* End if model includes a product */
+ else { /* no more sum */
+ /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
+ /* scanf("%d",i);*/
+ cutl(strd,strc,strb,'V');
+ ks++; /**< Number of simple covariates */
cptcovn++;
- cptcovprodnoage++;k1++;
- cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
- Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but
- because this model-covariate is a construction we invent a new column
- ncovcol + k1
- 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 */
- 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 */
- Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
- Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
- k2=k2+2;
- Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */
- Tvar[cptcovt+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */
- for (i=1; i<=lastobs;i++){
- /* Computes the new covariate which is a product of
- covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
- covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
- }
- } /* End age is not in the model */
- } /* End if model includes a product */
- else { /* no more sum */
- /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
- /* scanf("%d",i);*/
- cutl(strd,strc,strb,'V');
- ks++; /**< Number of simple covariates */
- cptcovn++;
- Tvar[k]=atoi(strd);
- }
- strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */
- /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
- scanf("%d",i);*/
- } /* end of loop + */
- } /* end model */
+ Tvar[k]=atoi(strd);
+ }
+ strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */
+ /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
+ scanf("%d",i);*/
+ } /* end of loop + on total covariates */
+ } /* end if strlen(modelsave == 0) age*age might exist */
+ } /* end if strlen(model == 0) */
/*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products.
If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/
}
ungetc(c,ficpar);
- fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
+ fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
numlinepar++;
- /* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */
- printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n",title, datafile, lastobs, firstpass,lastpass);
- /*
-
-
-
- */
- printf("\nftol=%e \n", ftol);
- printf("stepm=%d \n", stepm);
- printf("ncovcol=%d nlstate=%d \n", ncovcol, nlstate);
- printf("ndeath=%d maxwav=%d mle=%d weight=%d\n", ndeath, maxwav, mle, weightopt);
- printf("model=%s\n",model);
- fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
- fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+ printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
+ if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */
+ model[strlen(model)-1]='\0';
+ fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
+ fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
fflush(ficlog);
+ if(model[0]=='#'|| model[0]== '\0'){
+ printf("Error in 'model' line: model should start with 'model=1+age+' and end with '.' \n \
+ 'model=1+age+.' or 'model=1+age+V1.' or 'model=1+age+age*age+V1+V1*age.' or \n \
+ 'model=1+age+V1+V2.' or 'model=1+age+V1+V2+V1*V2.' etc. \n"); \
+ if(mle != -1){
+ printf("Fix the model line and run imach with mle=-1 to get a correct template of the parameter file.\n");
+ exit(1);
+ }
+ }
while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);
fgets(line, MAXLINE, ficpar);
v1+v2*age+v2*v3 makes cptcovn = 3
*/
if (strlen(model)>1)
- ncovmodel=2+nbocc(model,'+')+1; /*Number of variables including intercept and age = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/
+ ncovmodel=2+nbocc(model,'+')+1; /*Number of variables including intercept and age = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7,age*age makes 3*/
else
- ncovmodel=2;
- nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
+ ncovmodel=2; /* Constant and age */
nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */
npar= nforce*ncovmodel; /* Number of parameters like aij*/
if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){
/* Main decodemodel */
+
if(decodemodel(model, lastobs) == 1)
goto end;
nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);
ncodemax[1]=1;
Ndum =ivector(-1,NCOVMAX);
- if (ncovmodel > 2)
+ if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */
tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
/* Nbcode gives the value of the lth modality of jth covariate, in
V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/
}
/*--------- results files --------------*/
- fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
+ fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n");
git://henry.ined.fr / .git / commitdiff