--- imach/src/imach.c 2022/03/17 08:45:53 1.310
+++ imach/src/imach.c 2022/04/05 21:24:39 1.312
@@ -1,6 +1,18 @@
-/* $Id: imach.c,v 1.310 2022/03/17 08:45:53 brouard Exp $
+/* $Id: imach.c,v 1.312 2022/04/05 21:24:39 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ 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
@@ -1159,12 +1171,12 @@ typedef struct {
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.310 2022/03/17 08:45:53 brouard Exp $ */
+/* $Id: imach.c,v 1.312 2022/04/05 21:24:39 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
char version[]=__IMACH_VERSION__;
char copyright[]="March 2021,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, INED 2000-2021";
-char fullversion[]="$Revision: 1.310 $ $Date: 2022/03/17 08:45:53 $";
+char fullversion[]="$Revision: 1.312 $ $Date: 2022/04/05 21:24:39 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -3847,7 +3859,7 @@ double funcone( double *x)
/* Fixed */
/* 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<=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[1]]=covar[Tvar[TvarFind[1]]][i]; */
/* cov[2+6]=covar[Tvar[6]][i]; */
@@ -4702,10 +4714,13 @@ Title=%s
Datafile=%s Firstpass=%d La
if(s[m][iind]==-1)
printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
- for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean */
- idq[z1]=idq[z1]+weight[iind];
- 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]; /* *weight[iind];*/ /* Computes mean of quantitative with selected filter */
+ for (z1=1; z1<= nqfveff; z1++) { /* Quantitative variables, calculating mean on known values only */
+ if(!isnan(covar[ncovcol+z1][iind])){
+ idq[z1]=idq[z1]+weight[iind];
+ 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) */
/* printf("j=%d, j1=%d Age %d, iind=%d, num=%09ld m=%d\n",j,j1,(int)agev[m][iind],iind, num[iind],m); */
@@ -4768,16 +4783,19 @@ Title=%s
Datafile=%s Firstpass=%d La
Printing means of quantitative variables if any
*/
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]);
if(weightopt==1){
printf(" Weighted mean and standard deviation of");
fprintf(ficlog," 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]));
- 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]));
- fprintf(ficresphtmfr," 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}
+ var = \frac{\sum w (x-mu)^2}{\sum w} = \frac{w x^2}{\sum w} - mu^2
+ */
+ 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)
\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(m=1;m<=lastpass;m++){ */
@@ -5494,6 +5512,8 @@ void concatwav(int wav[], int **dh, int
if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */
switch(Fixed[k]) {
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*/
ij=(int)(covar[Tvar[k]][i]);
/* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
@@ -5507,8 +5527,10 @@ void concatwav(int wav[], int **dh, int
else if (ij < modmincovj)
modmincovj=ij;
if (ij <0 || ij >1 ){
- printf("Information, IMaCh doesn't treat covariate with missing values (-1), individual %d will be skipped.\n",i);
- fprintf(ficlog,"Information, currently 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,"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)){
printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
@@ -5583,6 +5605,16 @@ void concatwav(int wav[], int **dh, int
break;
} /* end switch */
} /* 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*/
for (k=-1; k< maxncov; k++) Ndum[k]=0;
@@ -9507,6 +9539,8 @@ int readdata(char datafile[], int firsto
cutv(stra, strb, line, ' ');
if(strb[0]=='.') { /* Missing value */
lval=-1;
+ coqvar[iv][i]=NAN;
+ covar[ncovcol+iv][i]=NAN; /* including qvar in standard covar for performance reasons */
}else{
errno=0;
/* what_kind_of_number(strb); */