All NAs produced when including all covariates, but runs fine when excluding one?

[Meghana Pagadala]

Hi!

I have been running plink for the past few months and have been developing my analyses to include some important covariates. I am currently working on running a logistic regression with a large sample of individuals (~10000) and only 175 cases. I have about 40 covariates, but only about 25 apply to the cases. So, I tried filtering the covariate file to the 25 covariates that apply. However, I keep getting the NAs (although this worked for other cases). I know this has something to do with the number of covariates/cases I have, but I am confused as to why the association runs fine when I exclude one of the covariates (have 24 total) but not with all of them. For 4 of the covariates, there is 1 case and several other controls. I am not sure if that is the reason. Any help would be appreciated!

Here is the command (plink v1.90b4, but does not work for plink2 either):

plink --bfile genotypes --chr 20 --covar all-tissue.cov --covar-name C1-C3, TISSUE1, TISSUE16, TISSUE20, TISSUE24, TISSUE17, TISSUE11, TISSUE18, TISSUE26, TISSUE1, TISSUE9, TISSUE28, TISSUE5, TISSUE25, TISSUE12, TISSUE7, TISSUE21, TISSUE15, TISSUE4, TISSUE30, TISSUE14, TISSUE0, TISUE23 --sex --logistic --adjust --pheno pheno.txt --out chr20

I'd be happy to send more info that would help to understand this better!

Thanks!

[Christopher Chang]

If TISSUE{0, 1, ..., 28} is a categorical covariate (i.e. all samples are a member of exactly one of those categories), you must remove one category (preferably the largest one) to avoid linear dependency between the categorical covariates and the intercept column.

[Meghana Pagadala]

Thank you for the quick response. 

So, would the best solution be then to exclude the categorical covariate with the largest number of cases or just overall number of samples?

Also, is there a reason why the logistic regression ran successfully for the other scenarios where I included all, but did not run successfully for this smaller case/control run?

Thanks again for the feedback!

[Christopher Chang]

Er, if the logistic regression worked in other scenarios with the same covariates, my initial statement about the categorical covariate wasn't quite right.

I recommend running the 31 Jul 2019 plink2 build with "--glm firth-fallback cols=+err"; this adds an extra output column which summarizes the reason behind each NA result.

[Meghana Pagadala]

Hi, Christopher.

I downloaded the 31 Jul build of plink2 and have been trying to run with the "--glm firth-fallback cols=+err" option, but it keeps giving me an error saying:

Error: Unrecognized ID 'err' in --glm column set descriptor.

I checked the build in the output and am sure I am running the correct version. Is there something that I am missing?

Thanks,

Meghana

[Christopher Chang]

Hmm, that's odd, I just re-downloaded the 31 Jul Linux AVX2 binary and didn't have a problem.  Maybe I screwed up one of the other binaries; what operating system/build are you on?

[Meghana Pagadala]

I'm on Linux 64-bit. I tried with the Linux 64-bit Intel buil and I could download plink2, just couldn't get the cols error to work.

Thanks,

Meghana

[Christopher Chang]

Hmm, that also works for me.  Can you post the .log file from a failed attempt to use cols=+err?

[Meghana Pagadala]

It worked this time. Perhaps I was calling something wrong.

I got the error message for why the P value is NA and it gives me this:

    This generates a fake input dataset with the specified number of samples

    and SNPs.

    * By default, the missing dosage and phenotype frequencies are zero.

      These can be changed by providing 3rd and 4th numeric parameters.

    * By default, allele codes are As and Bs; this can be changed with the

      'acgt', '1234', or '12' modifier.

    * By default, one binary phenotype is generated.  'pheno-ct=' can be used

      to change the number of phenotypes, and 'scalar-pheno' causes these

      phenotypes to be normally distributed scalars.

    * By default, all (nonmissing) dosages are in {0,1,2}.  To make some of

      them take on decimal values, use 'dosage-freq='.  (These dosages are

      affected by --hard-call-threshold and --dosage-erase-threshold.)

I guess could this be resolved by excluding the covariate with most cases?

Thanks,

Meghana

[Christopher Chang]

Er, can you clarify what was in the ERRCODE column?  It should be one of the following values:

  SAMPLE_CT<=PREDICTOR_CT

  CONST_ALLELE

  CORR_TOO_HIGH

  VIF_INFINITE

  VIF_TOO_HIGH

  SEPARATION

  RANK_DEFICIENT

  LOGISTIC_CONVERGE_FAIL

  FIRTH_CONVERGE_FAIL

  INVALID_RESULT

[Meghana Pagadala]

Oh, ya. 

It was LOGISTIC_CONVERGE_FAIL

Thanks,

Meghana

[Christopher Chang]

Is "firth-fallback" in the --glm part of your command line?

[Meghana Pagadala]

Oh, I think this one with the error didn't include it, because it was taking a really long time to run with the firth-fallback, but I have ones with the firth-fallback still running, but they are all at 0% writing status. I can wait until those are finished and see the file.

Thanks,

Meghana

[Meghana Pagadala]

With the firth-fallback, there seems to be no error. So, I guess that is the option to use?

Thanks,

Meghana

[Christopher Chang]

Yup.  This will become the default when alpha 3 builds are available.

[Meghana Pagadala]

Gotcha. So, I guess I can run this with all the covariates still?

Thanks,

Meghana

[Christopher Chang]

Yeah; plink2 will report an error when they actually aren't linearly independent.  You just seem to have separation, which Firth regression addresses.