Why does firth regression fail to converge?

Jonathan Margoliash

unread,

Feb 9, 2024, 12:02:35 AM2/9/24

to plink2-users

Hi there,

I'm running a GWAS on a case/control phenotype with 11k cases, 155k controls and 43 covariates . I'm trying to run with firth regression always on and cc-residualize and I get the error

Error: Cannot proceed with --glm regression on phenotype
'audit_combined_q4_binarized', since covariate-only Firth regression failed to
converge.

Any idea why this might happen? It doesn't seem like this is likely to be a power issue. Due to the size of the data, I'm intuiting I will need the performance boost this feature provides. I've posted the full stdout log below if it helps.

Thank you,

Jonathan Margoliash

--------

PLINK v2.00a6LM AVX2 Intel (5 Feb 2024) www.cog-genomics.org/plink/2.0/
(C) 2005-2024 Shaun Purcell, Christopher Chang GNU General Public License v3
Logging to plink2.log.
Options in effect:
--chr 4
--ci 0.99999995
--covar-name age is_from_followup sex pc1 pc2 pc3 pc4 pc5 pc6 pc7 pc8 pc9 pc10 pc11 pc12 pc13 pc14 pc15 pc16 pc17 pc18 pc19 pc20 pc21 pc22 pc23 pc24 pc25 pc26 pc27 pc28 pc29 pc30 pc31 pc32 pc33 pc34 pc35 pc36 pc37 pc38 pc39 pc40 --extract bed1 region.bed
--glm omit-ref pheno-ids hide-covar cols=+gcountcc,-tz,-nobs,-test single-prec-cc cc-residualize firth
--mac 20
--memory 56000
--no-psam-pheno
--pfile /cromwell-executions/gwass/d7f4486f-89cb-468b-80a7-9c093c8dad08/call-gwas/shard-0/gwas/2a882bca-9011-4c69-a735-5940777c2790/call-gwas/gwas/8b896998-9e30-4f8b-a26b-ec4e6f2a73cc/call-chromosomal_plink_snp_association/shard-82/inputs/-1087981777/chr4
--pheno /cromwell-executions/gwass/d7f4486f-89cb-468b-80a7-9c093c8dad08/call-gwas/shard-0/gwas/2a882bca-9011-4c69-a735-5940777c2790/call-gwas/gwas/8b896998-9e30-4f8b-a26b-ec4e6f2a73cc/call-chromosomal_plink_snp_association/shard-82/inputs/1805111232/out.tab
--pheno-name audit_combined_q4_binarized
--threads 28

Start time: Thu Feb 8 20:15:22 2024
257498 MiB RAM detected, ~237991 available; reserving 56000 MiB for main
workspace.
Using up to 28 threads (change this with --threads).
487409 samples (264296 females, 222987 males, 126 ambiguous; 487409 founders)
loaded from
/cromwell-executions/gwass/d7f4486f-89cb-468b-80a7-9c093c8dad08/call-gwas/shard-0/gwas/2a882bca-9011-4c69-a735-5940777c2790/call-gwas/gwas/8b896998-9e30-4f8b-a26b-ec4e6f2a73cc/call-chromosomal_plink_snp_association/shard-82/inputs/-1087981777/chr4.psam.
6555871 variants loaded from
/cromwell-executions/gwass/d7f4486f-89cb-468b-80a7-9c093c8dad08/call-gwas/shard-0/gwas/2a882bca-9011-4c69-a735-5940777c2790/call-gwas/gwas/8b896998-9e30-4f8b-a26b-ec4e6f2a73cc/call-chromosomal_plink_snp_association/shard-82/inputs/-1087981777/chr4.pvar.
1 binary phenotype loaded (10756 cases, 154663 controls).
--extract bed1: 6233317 variants excluded.
43 covariates loaded from /cromwell-executions/gwass/d7f4486f-89cb-468b-80a7-9c093c8dad08/call-gwas/shard-0/gwas/2a882bca-9011-4c69-a735-5940777c2790/call-gwas/gwas/8b896998-9e30-4f8b-a26b-ec4e6f2a73cc/call-chromosomal_plink_snp_association/shard-82/inputs/1805111232/out.tab.
Calculating allele frequencies... 0%^H^H17%^H^H^H37%^H^H^H58%^H^H^H78%^H^H^H98%^H^H^Hdone.
79786 variants removed due to allele frequency threshold(s)
(--maf/--max-maf/--mac/--max-mac).
242768 variants remaining after main filters.

End time: Thu Feb 8 20:15:34 2024

Christopher Chang

unread,

Feb 9, 2024, 12:11:47 AM2/9/24

to plink2-users

1. Does the problem remain if you remove the "single-prec-cc" modifier?

2. If yes, what happens if you try to perform Firth regression for that phenotype on your covariates in R?

Jonathan Margoliash

unread,

Feb 9, 2024, 1:24:04 PM2/9/24

to plink2-users

Thanks for the quick response.

1. single-prec-cc seems to be a requirement for running cc-residualize. Plink prints out the following:

Error: --glm 'cc-residualize'/'firth-residualize' requires 'hide-covar' and
'single-prec-cc' to be specified as well.

2. the logistf function from the R library logistf works once I recode the case control data from control=1, case=2 to control=0, case=1. I've printed the output below. I didn't time the function call, but it took more than 5 minutes and less than an hour, perhaps plink is erroring out because it takes too long?

> lf = logistf(data = data, formula = audit_combined_q4_binarized ~ age + is_from_followup + sex + pc1 + pc2 + pc3 + pc4 + pc5 + pc6 + pc7 + pc8 + pc9 + pc10 + pc11 + pc12 + pc13 + p[0/1850] + pc1 + pc2 + pc3 + pc4 + pc5 + pc6 + pc7 + pc8 + pc9 + pc10 + pc11 + pc12 + pc13 + pc14 + pc15 + pc16 + pc17 + pc18 + pc19 + pc20 + pc21 + pc22 + pc23 + pc24 + pc25 + pc26 + pc27 + pc28 + pc29 + pc30 + pc31 + pc32 + pc33 + pc34 + pc35 + pc36 + pc37 + pc38 + pc39 + pc40) > lf

logistf(formula = audit_combined_q4_binarized ~ age + is_from_followup +

sex + pc1 + pc2 + pc3 + pc4 + pc5 + pc6 + pc7 + pc8 + pc9 +

pc10 + pc11 + pc12 + pc13 + pc14 + pc15 + pc16 + pc17 + pc18 +

pc19 + pc20 + pc21 + pc22 + pc23 + pc24 + pc25 + pc26 + pc27 +

pc28 + pc29 + pc30 + pc31 + pc32 + pc33 + pc34 + pc35 + pc36 +

pc37 + pc38 + pc39 + pc40, data = data)

Model fitted by Penalized ML

Confidence intervals and p-values by Profile Likelihood

Coefficients:

(Intercept) age is_from_followup sex

-2.7699337185 -0.4646840600 0.0621125020 -0.1406630239

pc1 pc2 pc3 pc4

0.0067248589 -0.0131057745 -0.0070889247 0.0006665748

pc5 pc6 pc7 pc8

0.1107050310 0.0044602825 -0.0099517546 0.0146329649
pc9 pc10 pc11 pc12
-0.0116100187 0.0023048373 0.0744667945 -0.0209293447
pc13 pc14 pc15 pc16
0.0045605640 -0.0070380891 0.0001209848 0.0149019790
pc17 pc18 pc19 pc20
-0.0158339835 0.0101965738 -0.0028354683 -0.0051956460
pc21 pc22 pc23 pc24
-0.0104549635 -0.0034903288 0.0115261320 -0.0046322938
pc25 pc26 pc27 pc28
0.0014326025 0.0352398232 0.0028402224 0.0106282577
pc29 pc30 pc31 pc32
-0.0118609274 0.0189446382 -0.0045801335 -0.0001205553
pc33 pc34 pc35 pc36
-0.0049872798 0.0014844928 0.0064380209 0.0139228833
pc37 pc38 pc39 pc40
0.0044735929 0.0013223264 0.0095299207 -0.0020002599

Likelihood ratio test=2478.144 on 43 df, p=0, n=165419

Jonathan Margoliash

unread,

Feb 9, 2024, 1:37:24 PM2/9/24

to plink2-users

Two more pieces of information, which may or may not help:

plink2 firth regression still fails with the same error when I run it with the single covariate sex instead of all 43.

The error the logistf function in R's logistf package gives when I code the control/case data as 1/2 instead of 0/1 is

Error in logistf.fit(x = x, y = y, weight = weight, offset = offset, firth, :
In iteration 7: Determinant of Fisher information matrix was numerically 0

Christopher Chang

unread,

Feb 9, 2024, 2:05:39 PM2/9/24

to plink2-users

Oops, I forgot about the "cc-residualize requires single-prec-cc" detail, my bad.

Is it possible for you to post or send me just the phenotype and sex columns that are sufficient to trigger the failure?

On Friday, February 9, 2024 at 10:24:04 AM UTC-8 Jonathan Margoliash wrote:

Jonathan Margoliash

unread,

Feb 9, 2024, 2:16:39 PM2/9/24

to plink2-users

Attached. I've removed the IID and FID columns as well as the other covariates so as to not expose information. Let me know if there's anything else I can provide.

test_firth_pheno_public.tab

Christopher Chang

unread,

Feb 9, 2024, 5:19:31 PM2/9/24

to plink2-users

Thanks. Looks like with that data, single precision math simply isn't good enough for the Firth regression loop to meet the usual convergence criteria, it gets stuck in a loop revolving around the solution.

I'm not inclined to second-guess the convergence criteria, so I will try to implement double-precision cc-residualize later this month.

Jonathan Margoliash

unread,

Feb 9, 2024, 9:32:00 PM2/9/24

to plink2-users

Okay, good to know why that's happening. Thank you for all the prompt assistance!

Message has been deleted

Chris Chang

unread,

Mar 31, 2024, 11:55:20 PM3/31/24

to Jonathan Margoliash, plink2-users

fyi, cc-residualize is usable without single-prec-cc in the 18 Mar build.

--
You received this message because you are subscribed to the Google Groups "plink2-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plink2-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plink2-users/1fe125aa-ea3a-42c6-ba27-69cabf771e74n%40googlegroups.com.

David Jakubosky

unread,

Sep 17, 2024, 8:18:52 PM9/17/24

to plink2-users

I'm not sure if this is the right place to reply - but I am still unable to get cc-residualize working without the single-prec-cc flag in the most recent build (20 Aug 2024) - am I missing something obvious?

PLINK v2.00a5.14LM 64-bit Intel (20 Aug 2024) www.cog-genomics.org/plink/2.0/

Logging to /local_disk0/chr22_test.log.
Options in effect:
--bfile /local_disk0/chr22_variants
--covar-name EstimatedAge,Sex,PC1,PC2,PC3,PC4,PC5
--covar-variance-standardize
--glm firth-fallback hide-covar cc-residualize
--mac 25
--no-input-missing-phenotype
--out /local_disk0/chr22_test
--pheno /local_disk0/phenotypes
--pheno-name pheno_name.txt

Start time: Wed Sep 18 00:00:55 2024

Error: --glm 'cc-residualize'/'firth-residualize' requires 'hide-covar' and
'single-prec-cc' to be specified as well.

Christopher Chang

unread,

Sep 17, 2024, 8:20:33 PM9/17/24

to plink2-users

Yes, that wasn't a trivial code change so you need to use the development build (a6), not a5.

Reply all

Reply to author

Forward