variani / lme4qtl

Mixed models @lme4 + custom covariances + parameter constraints
GNU General Public License v3.0
48 stars 10 forks source link

PIRLS step-halvings failed to reduce deviance in pwrssUpdate with probit relmatGlmer #30

Closed zoeluo15 closed 2 years ago

zoeluo15 commented 2 years ago

Hi, I am using the relmatGlmer function with the probit link function in lme4qtl. Here is my data and code

kinship.txt genotype_phenotype.txt

test<-read.table("genotype_phenotype.txt")
kinship<-as.matrix(read.table("kinship.txt")) colnames(kinship)<-rownames(kinship) relmatGlmer(phenotype ~ genotype + (1|ID), test, relmat = list(ID = kinship), family = binomial(link="probit"))

And I got the following error message:

Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev, : (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

I found a relevant question in lme4 https://github.com/lme4/lme4/issues/579, and tried their solution

relmatGlmer(phenotype ~ genotype + (1|ID), test, relmat = list(ID = kinship), family = binomial(link="probit"),nAGQ=20)

It returns the error message:

Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev, : AGQ only defined for a single scalar random-effects term

I think my random effect term is a scalar random effect. Am I doing anything wrong here?

Any help would be appreciated!

Zoe

variani commented 2 years ago

Hi Zoe,

I'm not sure what is exactly happening here. Looking at the issue for lme4 (you mentioned), the cause seems to be in a low number of observations. I wonder if you can simulate data, increase the number of observations and check if the issue is gone.

Best, Andrey

zoeluo15 commented 2 years ago

Hi Andrey,

There is indeed no problem once I increase the number of observations in each group. But I would like to fit the model with an individual random effect (i.e. one observation per level of the random effect) because this is actually a real bird dataset.

I think it might be the distribution of the genotype (the fixed effect) as there is no problem fitting the same model with genotypes at another position (test2.txt).

table(test$genotype) #this is the dataset I attached in the original post

0 1 2 82 64 8

table(test2$genotype) #this is the new one that has no problem in model fitting

0 1 2 136 17 1

Many thanks, Zoe

variani commented 2 years ago

Hi Zoe,

Unfortunately, I cannot help too much here, as I was not the person who implemented the model fitting algorithm involved.

Perhaps, you can play with glmerControl. In particular, there are two algorithms for the optimizer (the first argument of glmerControl).

Best, Andrey

zoeluo15 commented 2 years ago

Thanks! I will have a try with glmerControl. - Zoe