melff / mclogit Goto Github PK

View Code? Open in Web Editor NEW

22.0 22.0 4.0 4.38 MB

mclogit: Multinomial Logit Models, with or without Random Effects or Overdispersion

Home Page: http://melff.github.io/mclogit/

R 100.00%

mclogit's People

Contributors

Stargazers

Watchers

Forkers

sfahnens pmcharrison yalchik rvlenth

mclogit's Issues

Syntax for increasing iterations of mblogit model

What does the control syntax look like for increasing the number of iterations of an mblogit (or mclogit) model?

I have tried:

lm.pred.da.RE.subj <- mblogit(thesis_da ~ log(time_since_last_utt), 
                                  random = c(~1|subj_id),
                                  control= mclogit.control(list(maxit=50,epsilon=.001)),
                             data=study1.dat.filter)

But I get the error:
Error in mclogit.control(list(maxit = 50, epsilon = 0.001)) : value of epsilon must be > 0

I am attempting to see if my model will converge after a few more iterations. I know it may just be an issue with my model, though.

rename residual.df to df.residual in results object

This is only a very minor thing but the results object from a mclogit call should return an object where the residual df is called "df.residual". At the moment (0.6.1) it is called "residual.df" which can produce errors if called by another function that uses the normal regression result output object structure.

[Question] Interpreting and plotting a mblogit model

Dear Professor Elff,

First of all, thank you very much for the package you published which is a perfect tool to analyze my data! I am currently working on my thesis concerning emotional speech in different languages. Specifically, the goal of the study is to confirm which acoustic parameter is significantly related to each emotion. Very below of this post are the models I built solar for English and Korean data.

My questions are as follows:

Since the reference level for the English model is 'ang', 'hap' and 'sad' were respectively compared as shown in the result. Is there a way to know the relationship between 'hap' and 'sad' in this case? When I built another model having 'hap' as a reference level, exp(coef(model)) was different from that of the model where 'ang' was assigned as a reference level.
(Co)variance: Is it normal to have a large (co)variance value? Unlike the English model, the Korean model showed high numbers of (Co)variance as below (The model is not converged). Does it imply that the random effects are innegligible? What do the numbers exactly mean?
Is it possible to plot the model reflecting the mixed effect? If so, could you please provide me some guidance for it? (I have tried with plot(eng_m2) and it did not work...)
How can I know which model fits the best? I saw at a book that anova is not recommended for this package, and you also mentioned somewhere here that bic would not be good either to trust.
Lastly, as I am relatively a beginner in statistics, I am not sure whether I have to check multicollinearity for my models. When I tries, R gives warning message when I entered 'vif(eng_m2) as 'Warning message:In vif.default(eng_m2) : No intercept: vifs may not be sensible.'.

Thank you so much for reading the questions!
Best regards,
Jueun Kang

English model

> str(eng_scaled)
tibble [18,231 × 11] (S3: tbl_df/tbl/data.frame)
 $ lang     : Factor w/ 1 level "eng": 1 1 1 1 1 1 1 1 1 1 ...
 $ gender   : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 1 1 1 1 ...
 $ emotion  : Factor w/ 3 levels "ang","hap","sad": 1 1 1 1 1 1 1 1 1 1 ...
 $ speaker  : Factor w/ 10 levels "eng1","eng10",..: 7 7 7 7 7 7 7 7 7 7 ...
 $ phone    : Factor w/ 3 levels "a","i","u": 3 3 3 3 3 3 3 3 3 3 ...
 $ F0       : num [1:18231] 1.22 0.88 1.71 0.49 0.65 1.05 1.92 1.44 1.64 0.68 ...
 $ F1       : num [1:18231] 0.45 0.23 -0.82 -0.01 -1.02 -1.08 -0.82 -0.71 -0.73 0.1 ...
 $ F2       : num [1:18231] -1.15 -1.51 -1.53 -1.56 0.62 -1.87 -1.78 -0.19 -1.55 -1.41 ...
 $ duration : num [1:18231] -0.16 -0.43 -0.43 -0.84 -1.25 -0.43 -0.16 -1.39 -1.12 -0.98 ...
 $ intensity: num [1:18231] 1.22 0.02 0.62 -0.58 -0.87 0.62 2.56 1.07 1.37 1.52 ...
 $ energy   : num [1:18231] 0.4 -0.23 -0.29 -0.55 -0.61 0.01 5.19 -0.39 -0.03 0.25 ...

> summary(eng_m2)
Call:
mblogit(formula = emotion ~ F0 + F1 + F2 + duration + intensity + 
    energy, data = eng_scaled, random = list(~1 | speaker, ~1 | 
    gender, ~1 | phone))

Equation for hap vs ang:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.15282    0.33990  -0.450    0.653    
F0           0.76615    0.02707  28.297  < 2e-16 ***
F1          -0.25688    0.03020  -8.506  < 2e-16 ***
F2          -0.16723    0.03202  -5.222 1.77e-07 ***
duration     0.27422    0.02360  11.620  < 2e-16 ***
intensity   -0.17466    0.02942  -5.937 2.91e-09 ***
energy      -0.58378    0.03440 -16.971  < 2e-16 ***

Equation for sad vs ang:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.20731    0.44908  -0.462   0.6443    
F0          -1.09574    0.03727 -29.402  < 2e-16 ***
F1           0.05779    0.02996   1.929   0.0538 .  
F2           0.22005    0.03359   6.551 5.71e-11 ***
duration     0.21687    0.02359   9.192  < 2e-16 ***
intensity    0.25189    0.03087   8.160 3.36e-16 ***
energy      -0.64147    0.04042 -15.871  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Co-)Variances:
Grouping level: speaker 
      Estimate            Std.Err.         
hap~1  0.14634            0.003384         
sad~1 -0.09344  0.20440   0.004423 0.006064

>Grouping level: gender 
      Estimate            Std.Err.         
hap~1  0.16080            0.007254         
sad~1 -0.07084  0.27069   0.011991 0.023856

Grouping level: phone 
      Estimate              Std.Err.           
hap~1 0.0594145             0.0001714          
sad~1 0.0002151 0.1353236   0.0004175 0.0020331

Null Deviance:     40060 
Residual Deviance: 35200 
Number of Fisher Scoring iterations:  6
Number of observations
  Groups by speaker: 10
  Groups by gender: 2
  Groups by phone: 3
  Individual observations:  18231

Korean model

> summary(kor_m1)

Call:
mblogit(formula = emotion ~ F0 + F1 + F2 + duration + intensity + 
    energy, data = kor_scaled, random = list(~1 | speaker, ~1 | 
    gender, ~1 | phone))

Equation for hap vs ang:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.69225    4.00926   0.422  0.67296    
F0           0.77587    0.13249   5.856 4.74e-09 ***
F1          -0.66141    0.10112  -6.541 6.12e-11 ***
F2          -0.65357    0.15319  -4.266 1.99e-05 ***
duration     0.15896    0.09844   1.615  0.10637    
intensity   -0.26378    0.09987  -2.641  0.00826 ** 
energy      -0.25351    0.13002  -1.950  0.05121 .  

Equation for sad vs ang:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.0763     6.3832  -0.169 0.866098    
F0           -0.4551     0.2555  -1.781 0.074872 .  
F1           -0.9775     0.3521  -2.776 0.005498 ** 
F2           -0.2067     0.2924  -0.707 0.479654    
duration      0.1314     0.2177   0.604 0.546162    
intensity    -0.9113     0.2696  -3.380 0.000725 ***
energy        0.3450     0.2803   1.231 0.218363    
---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Co-)Variances:
Grouping level: speaker 
      Estimate        Std.Err.   
hap~1 105.71          142.5      
sad~1  90.06 222.87   195.8 370.3

Grouping level: gender 
      Estimate        Std.Err.   
hap~1  7.118          16.77      
sad~1 -6.603  6.818   16.44 16.06

Grouping level: phone 
      Estimate      Std.Err.   
hap~1 18.73          61.3      
sad~1 28.08 74.67   124.7 288.6

Null Deviance:     31510 
Residual Deviance: 1594 
Number of Fisher Scoring iterations:  25
Number of observations
  Groups by speaker: 20
  Groups by gender: 2
  Groups by phone: 3
  Individual observations:  14339
Note: Algorithm did not converge.

Problems with model definition

Thank you, the library has been very useful to me.

My problem is when I want to fit a model from the definition of formula (as.formula ()) and also add random effects:

# this worked
m <-mclogit::mblogit(y ~ x1+x2,random = ~1|id,
                 data = df)

# this did not work
model<-as.formula(y~x1+x2)
m <-mclogit::mblogit(model,random = ~1|id,
                 data = df)

output: Error in attributes(.Data) <- c(attributes(.Data), attrib) : 
  cannot set attribute on a symbol

The reason to define the model like this is that I am creating a function that calls to mblogit.

Thank you for your time.

systematic error with nested model

Hi, I really like this package as I find it intuitive and great for multinomial regression - thanks for building and sharing it. I am having an issue with using it for nested designs though. I have tried the following implementation (illustrated here with a toy dataset) on a number of datasets and I'm systematically getting a "Error in *tmp*[[k]] : subscript out of bounds" error (right after the statement "converged" gets printed in the trace). Is there something I'm missing, or a workaround perhaps? Best wishes, E.

set.seed(1)
country = c("france","ireland")
county = c("cork","dublin","limerick","galway")
dept = c("var","alpes-maritimes","vaucluse","bouches-du-rhone")
fake.fr = expand.grid(country[1],dept)
fake.ie = expand.grid(country[2],county)
for(i in 1:6){ # creating 512 rows of this data
fake.fr = rbind(fake.fr,fake.fr)
fake.ie = rbind(fake.ie,fake.ie)
}
fake = rbind(fake.fr,fake.ie)
n = nrow(fake)
fake$x1 = runif(n)
fake$x2 = rnorm(n)
fake$y = factor(c("coffee","tea")[sample(1:2,size=n,replace=T)])
names(fake)[1:2] = c("country","county")

level-2 random effects works:

(om = mblogit(y~x1+x2, random=~1|country, data=fake))

level-2 and level-3 random effects does not work:

(om = mblogit(y~x1+x2, random=~1|county/country, data=fake))

can I conduct a likelihood ratio test with output from a mixed effects model in mblogit?

Hi.

I have been fitting fixed effects models with mblogit and then conducting likelihood ratio tests using the anova command and the pchisq command.

Today, my first mixed effects model converged and I checked to see which method is better to use for conducting a likelihood ratio test: MQL or PQL. Unfortunately, the information I read suggests neither.

Can you please tell me if there is a way to conduct a likelihood ratio test with the output from a mixed effects model in mblogit, and if not, can you suggest a test that will work. I am testing the effect of a single treatment, in the presence of a random effect of individual, and sometimes with other fixed covariates. My response variable has three categories, and I am testing for an effect in both logits.

I will paste some sample lines of code below.

Thank you very much!

Rebecca

MBLrWM <- mblogit(aBehavior ~ Trt15_20_10min, random=~1|Deployment, data=Swm, maxit = 100, method="MQL")
MBLrW0M <- mblogit(aBehavior ~ 1, random=~1|Deployment, data=Swm, maxit = 100, method="MQL")
anova(MBLrW0M,MBLrWM)
pchisq(0.53241,df=2,lower.tail = FALSE)

Frequent error in mblogit "Error in solve.default(Information) : system is computationally singular: reciprocal condition number" with models that run OK in nnet::multinom

I am frequently encountering fitting errors in mblogit, returning me the error
"Error in solve.default(Information) : system is computationally singular: reciprocal condition number"
for models that run OK in nnet::multinom. Would you happen to have any idea on how to resolve these?
A reproducible example is given below. It could be that there is some problem with complete separation / some Hauck-Donner effect, but given that it happens so frequently, also in situations where the same fit in nnet::multinom runs OK I have the feeling there might be some bug in mblogit. Either way, any advice on how to get rid of this error would be welcome! In logistic regression I recall that one can add a small amount of noise to zero counts or add some ridge regularisation (e.g. by row augmenting the covariate matrix with a diagonal matrix with sqrt(lambda) with lambda=ridge penalty) to get rid of this sort of problems - I don't know if this could be a solution also for multinomial regression?

# data=SARS-CoV2 coronavirus variants (variant) through time (collection_date_num)
# in India, count=actual count (nr of sequenced genomes)

dat = read.csv("https://www.dropbox.com/s/u27cn44p5srievq/dat.csv?dl=1")
dat$variant = factor(dat$variant)

# 1. multinom::net MULTINOMIAL FIT ####
library(nnet)
library(splines)
set.seed(1)
fit_nnet = nnet::multinom(variant ~ ns(collection_date_num, df=2), 
                          weights=count, data=dat)
summary(fit_nnet)
# Coefficients:
#   (Intercept) ns(collection_date_num, df = 2)1
# Beta                 3.798225                        4.9204017
# Delta              -10.043973                       49.4947796
# Omicron (BA.1)    -265.410841                      487.1385690
# Omicron (BA.2)    -467.537309                      831.1140917
# Omicron (BA.2.74)   -4.616715                       15.4357149
# Omicron (BA.2.75)   11.365092                      -38.7658668
# Omicron (BA.2.76)    5.437985                        0.3403524
# Omicron (BA.3)     -15.924342                       42.5896877
# Omicron (BA.4)      -7.440555                       23.8097760
# Omicron (BA.5)      -2.940342                       15.3461333
# Other               19.322779                       -8.4319599
# ns(collection_date_num, df = 2)2
# Beta                                      35.17926
# Delta                                     55.21257
# Omicron (BA.1)                           195.27695
# Omicron (BA.2)                           312.90170
# Omicron (BA.2.74)                         75.53290
# Omicron (BA.2.75)                         79.80443
# Omicron (BA.2.76)                         70.73509
# Omicron (BA.3)                            75.60153
# Omicron (BA.4)                            73.26847
# Omicron (BA.5)                            74.45010
# Other                                     54.14069
# 
# Std. Errors:
#   (Intercept) ns(collection_date_num, df = 2)1
# Beta                1.8548783                        1.9750676
# Delta               0.8914081                        0.7919794
# Omicron (BA.1)      8.9955175                       15.6207548
# Omicron (BA.2)      7.0504753                       12.1069603
# Omicron (BA.2.74)   1.1513095                        0.7721499
# Omicron (BA.2.75)   2.5900436                        5.7317343
# Omicron (BA.2.76)  15.9802663                       27.2148132
# Omicron (BA.3)     27.1955306                       46.8288796
# Omicron (BA.4)      1.0868226                        0.6370319
# Omicron (BA.5)      1.3766617                        1.5280076
# Other               0.8118766                        0.6288743
# ns(collection_date_num, df = 2)2
# Beta                                      5.998582
# Delta                                     3.133279
# Omicron (BA.1)                            5.388977
# Omicron (BA.2)                            4.888832
# Omicron (BA.2.74)                         3.177189
# Omicron (BA.2.75)                         3.862689
# Omicron (BA.2.76)                         9.112225
# Omicron (BA.3)                           14.632091
# Omicron (BA.4)                            3.116746
# Omicron (BA.5)                            3.081921
# Other                                     3.128327
# 
# Residual Deviance: 128392.2 
# AIC: 128458.2 

# predicted proportions at last day calculated using emmeans:
library(emmeans)
multinom_emmeans = emmeans(fit_nnet, ~ variant,  
                       mode = "prob",
                       at=list(collection_date_num = 
                                 max(data_agbyweek1$collection_date_num)))
multinom_emmeans
# variant               prob       SE df lower.CL upper.CL
# Alpha             0.00e+00 0.00e+00 33 0.00e+00 0.00e+00
# Beta              0.00e+00 0.00e+00 33 0.00e+00 0.00e+00
# Delta             7.73e-06 1.17e-06 33 5.34e-06 1.01e-05
# Omicron (BA.1)    1.82e-04 6.42e-05 33 5.14e-05 3.13e-04
# Omicron (BA.2)    1.76e-01 7.45e-03 33 1.61e-01 1.91e-01
# Omicron (BA.2.74) 9.03e-02 7.98e-03 33 7.41e-02 1.07e-01
# Omicron (BA.2.75) 1.68e-01 1.90e-02 33 1.30e-01 2.07e-01
# Omicron (BA.2.76) 2.89e-01 1.35e-02 33 2.62e-01 3.16e-01
# Omicron (BA.3)    1.34e-02 2.10e-03 33 9.10e-03 1.76e-02
# Omicron (BA.4)    1.67e-02 2.47e-03 33 1.17e-02 2.17e-02
# Omicron (BA.5)    2.03e-01 1.08e-02 33 1.81e-01 2.25e-01
# Other             4.23e-02 3.15e-03 33 3.59e-02 4.87e-02
# 
# Confidence level used: 0.95 

# 2. SAME MULTINOMIAL FIT USING mblogit ####
library(mclogit)
fit_mblogit = mblogit(variant ~ ns(collection_date_num, df=2),
                      weight=count,
                      data=dat,
                      from.table=TRUE, dispersion=FALSE,
                      control=mclogit.control(maxit=27))
summary(fit_mblogit)
# Equation for Beta vs Alpha:
#   Estimate Std. Error z value Pr(>|z|)
# (Intercept)                      -2.140e+01  1.505e+06       0        1
# ns(collection_date_num, df = 2)1 -9.047e+00  3.062e+06       0        1
# ns(collection_date_num, df = 2)2  2.608e+01  1.379e+06       0        1
# 
# Equation for Delta vs Alpha:
#   Estimate Std. Error z value Pr(>|z|)
# (Intercept)                      -2.140e+01  1.505e+06       0        1
# ns(collection_date_num, df = 2)1 -9.047e+00  3.062e+06       0        1
# ns(collection_date_num, df = 2)2  2.608e+01  1.379e+06       0        1
# 
# Equation for Omicron (BA.1) vs Alpha:
#   Estimate Std. Error z value
# Omicron (BA.1)~(Intercept)                      -2.140e+01  1.505e+06       0
# Omicron (BA.1)~ns(collection_date_num, df = 2)1 -9.047e+00  3.062e+06       0
# Omicron (BA.1)~ns(collection_date_num, df = 2)2  2.608e+01  1.379e+06       0
# Pr(>|z|)
# Omicron (BA.1)~(Intercept)                             1
# Omicron (BA.1)~ns(collection_date_num, df = 2)1        1
# Omicron (BA.1)~ns(collection_date_num, df = 2)2        1
# 
# Equation for Omicron (BA.2) vs Alpha:
#   Estimate Std. Error z value
# Omicron (BA.2)~(Intercept)                      -6.051e+00  5.785e+05       0
# Omicron (BA.2)~ns(collection_date_num, df = 2)1 -3.120e+01  1.257e+06       0
# Omicron (BA.2)~ns(collection_date_num, df = 2)2  5.841e+01  3.135e+05       0
# Pr(>|z|)
# Omicron (BA.2)~(Intercept)                             1
# Omicron (BA.2)~ns(collection_date_num, df = 2)1        1
# Omicron (BA.2)~ns(collection_date_num, df = 2)2        1
# 
# Equation for Omicron (BA.2.74) vs Alpha:
#   Estimate Std. Error
# Omicron (BA.2.74)~(Intercept)                      -7.065e+00  5.791e+05
# Omicron (BA.2.74)~ns(collection_date_num, df = 2)1 -2.952e+01  1.258e+06
# Omicron (BA.2.74)~ns(collection_date_num, df = 2)2  5.573e+01  3.116e+05
# z value Pr(>|z|)
# Omicron (BA.2.74)~(Intercept)                            0        1
# Omicron (BA.2.74)~ns(collection_date_num, df = 2)1       0        1
# Omicron (BA.2.74)~ns(collection_date_num, df = 2)2       0        1
# 
# Equation for Omicron (BA.2.75) vs Alpha:
#   Estimate Std. Error
# Omicron (BA.2.75)~(Intercept)                      -6.311e+00  5.785e+05
# Omicron (BA.2.75)~ns(collection_date_num, df = 2)1 -3.077e+01  1.257e+06
# Omicron (BA.2.75)~ns(collection_date_num, df = 2)2  5.771e+01  3.131e+05
# z value Pr(>|z|)
# Omicron (BA.2.75)~(Intercept)                            0        1
# Omicron (BA.2.75)~ns(collection_date_num, df = 2)1       0        1
# Omicron (BA.2.75)~ns(collection_date_num, df = 2)2       0        1
# 
# Equation for Omicron (BA.2.76) vs Alpha:
#   Estimate Std. Error
# Omicron (BA.2.76)~(Intercept)                      -6.072e+00  5.784e+05
# Omicron (BA.2.76)~ns(collection_date_num, df = 2)1 -3.116e+01  1.257e+06
# Omicron (BA.2.76)~ns(collection_date_num, df = 2)2  5.835e+01  3.137e+05
# z value Pr(>|z|)
# Omicron (BA.2.76)~(Intercept)                            0        1
# Omicron (BA.2.76)~ns(collection_date_num, df = 2)1       0        1
# Omicron (BA.2.76)~ns(collection_date_num, df = 2)2       0        1
# 
# Equation for Omicron (BA.3) vs Alpha:
#   Estimate Std. Error z value
# Omicron (BA.3)~(Intercept)                      -2.140e+01  1.505e+06       0
# Omicron (BA.3)~ns(collection_date_num, df = 2)1 -9.047e+00  3.062e+06       0
# Omicron (BA.3)~ns(collection_date_num, df = 2)2  2.608e+01  1.379e+06       0
# Pr(>|z|)
# Omicron (BA.3)~(Intercept)                             1
# Omicron (BA.3)~ns(collection_date_num, df = 2)1        1
# Omicron (BA.3)~ns(collection_date_num, df = 2)2        1
# 
# Equation for Omicron (BA.4) vs Alpha:
#   Estimate Std. Error z value
# Omicron (BA.4)~(Intercept)                      -7.442e+00  5.800e+05       0
# Omicron (BA.4)~ns(collection_date_num, df = 2)1 -2.890e+01  1.260e+06       0
# Omicron (BA.4)~ns(collection_date_num, df = 2)2  5.475e+01  3.110e+05       0
# Pr(>|z|)
# Omicron (BA.4)~(Intercept)                             1
# Omicron (BA.4)~ns(collection_date_num, df = 2)1        1
# Omicron (BA.4)~ns(collection_date_num, df = 2)2        1
# 
# Equation for Omicron (BA.5) vs Alpha:
#   Estimate Std. Error z value
# Omicron (BA.5)~(Intercept)                      -6.825e+00  5.785e+05       0
# Omicron (BA.5)~ns(collection_date_num, df = 2)1 -2.992e+01  1.257e+06       0
# Omicron (BA.5)~ns(collection_date_num, df = 2)2  5.635e+01  3.120e+05       0
# Pr(>|z|)
# Omicron (BA.5)~(Intercept)                             1
# Omicron (BA.5)~ns(collection_date_num, df = 2)1        1
# Omicron (BA.5)~ns(collection_date_num, df = 2)2        1
# 
# Equation for Other vs Alpha:
#   Estimate Std. Error z value Pr(>|z|)
# (Intercept)                      -2.140e+01  1.505e+06       0        1
# ns(collection_date_num, df = 2)1 -9.047e+00  3.062e+06       0        1
# ns(collection_date_num, df = 2)2  2.608e+01  1.379e+06       0        1
# 
# Null Deviance:     13540 
# Residual Deviance: 7.354e-10 
# Number of Fisher Scoring iterations:  27 
# Number of observations:  43160 
# 
# 
# Note: Algorithm did not converge.

As you can see this fit returns unusually large standard errors, and increasing maxit to >30 results in the error

# Error in solve.default(Information) : 
#  system is computationally singular: reciprocal condition number = 1.66419e-16

It could be that this is due to complete separation / some Hauck-Donner effect, but given that this model runs OK in nnet::multinom I am tempted to think it could rather be a bug in mblogit.

# predicted variant proportions on last day calculated using emmeans :
mblogit_emmeans = emmeans(fit_mblogit, ~ variant,  
                           mode = "prob",
                           at=list(collection_date_num = 
                                     max(data_agbyweek1$collection_date_num)))
mblogit_emmeans
# variant             prob       SE  df asymp.LCL asymp.UCL
# Alpha             0.0000 9.00e-08 Inf -1.80e-07  2.00e-07
# Beta              0.0000 1.00e-08 Inf -2.00e-08  0.00e+00
# Delta             0.0000 1.00e-08 Inf -2.00e-08  0.00e+00
# Omicron (BA.1)    0.0000 1.00e-08 Inf -2.00e-08  0.00e+00
# Omicron (BA.2)    0.3653 3.73e-02 Inf  2.92e-01  4.38e-01
# Omicron (BA.2.74) 0.0299 1.32e-02 Inf  4.09e-03  5.58e-02
# Omicron (BA.2.75) 0.1916 3.05e-02 Inf  1.32e-01  2.51e-01
# Omicron (BA.2.76) 0.3473 3.68e-02 Inf  2.75e-01  4.20e-01
# Omicron (BA.3)    0.0000 1.00e-08 Inf -2.00e-08  0.00e+00
# Omicron (BA.4)    0.0120 8.42e-03 Inf -4.52e-03  2.85e-02
# Omicron (BA.5)    0.0539 1.75e-02 Inf  1.96e-02  8.81e-02
# Other             0.0000 1.00e-08 Inf -2.00e-08  0.00e+00
# 
# Confidence level used: 0.95

# NOTE: does not quite match emmeans nnet::multinom output
# above, but maybe because I couldn't run mblogit until convergence

Working covariance matrix

Hi Prof. Martin Elff,
I looked into your code of function mmclogit.fitPQLMQL, and found that you define the working covariance matrix W the same as the variance and covariance of multinomial distribution. W <- Diagonal(x=w*pi)-tcrossprod(W). I am curious about the reason why. Because according to the paper of Breslow and Clayton in 1993, W is defined as $(\phia_iv(\mu_i)*g'(\mu_i)^2)^{-1}$, where v(\mu_i) is 'W' in your code. It appears that now, your W^{-1} is not the variance covariance matrix of the working response variable y.star.

Thanks in advance.

Error "Error: no valid set of coefficients has been found: please supply starting values" for model that runs OK in nnet::multinom

Recently bumped into a model for which no valid set of coefficients could be found, even though the model runs OK in nnet::multinom :

library(nnet)
library(mclogit)
library(splines)
data(mtcars)
dat = mtcars
dat$cyl = factor(dat$cyl)
dat$am = factor(dat$am)

fit_multinom = multinom(cyl ~ ns(mpg, df=2) + am, data = dat) # runs OK

fit_mblogit = mblogit(cyl ~ ns(mpg, df=2) + am, data = dat, 
                       from.table=FALSE) 
# returns error "Error: no valid set of coefficients has been found: please supply starting values"

Is there any way to pass starting coefficients actually?

Including multiple random effects in an mblogit model

Hello,

I am wondering if it is possible to fit multiple un-nested random effects within an mblogit model. I am able to run the model with either random effect and can nest them as well, but in this case it does not make sense to nest the random effects and would be best to treat them separately. For example, in lmer I know the syntax would be y ~ x + (1|random effect 1) + (1|random effect 2), but is it possible to do something similar to this in an mblogit model? Happy to provide additional information and code but thought it would be best to ask this question at a conceptual level first before delving into specifics.

Thank you!

Algorithm did not converge

Hi Prof. Martin Elff,

I am trying to use your mblogit function to fit a multinomial logistic regression to examine the factors which predict participants in my experiment selecting one of three answers. Because my model is mixed (I am comparing different age groups of participants and also participants' performance on one of two trial types) I am including subject ID as a random factor.

For model selection, my approach is to start by fitting a full model (3 IVs with all their interactions) then removing components and comparing the BIC score of the model before and after removing each component. If removing a component dramatically reduces the BIC score (based on Raftery, 1995) then I will argue it does not add predictive value to the model and that component will not be included in the final model.

My problem here is that with the full model R studio returns the error "Algorithm does not converge". I increased the number of iterations to 200 but it returns the following issue:

> study1_fullModel <- mblogit(Answer ~ Age_group:Relief:Story, 
                                data = long_data2,
                                maxit = 200,
                                method = "PQL",
                                trace = T,
                                random = ~1|ID)
Iteration 1 - deviance = 866.3455 - criterion = 0.6663194
Iteration 2 - deviance = 666.6792 - criterion = 0.07784625
Iteration 3 - deviance = 538.5546 - criterion = 0.04248887
Iteration 4 - deviance = 562.2433 - criterion = 0.02415636
Iteration 5 - deviance = 584.9369 - criterion = 0.02079293
Iteration 6 - deviance = 592.4505 - criterion = 0.0183108
Iteration 7 - deviance = 585.5645 - criterion = 0.01538235
Iteration 8 - deviance = 569.38 - criterion = 0.01273065
Iteration 9 - deviance = 548.7899 - criterion = 0.01054194
Iteration 10 - deviance = 527.194 - criterion = 0.008786208
Iteration 11 - deviance = 506.3881 - criterion = 0.007371137
Iteration 12 - deviance = 487.2235 - criterion = 0.006213371
Iteration 13 - deviance = 470.1385 - criterion = 0.005256752
Iteration 14 - deviance = 455.3552 - criterion = 0.004465433
Iteration 15 - deviance = 442.9199 - criterion = 0.003814404
Iteration 16 - deviance = 432.7273 - criterion = 0.003283286
Iteration 17 - deviance = 424.5586 - criterion = 0.002853174
Iteration 18 - deviance = 418.1291 - criterion = 0.002505815
Iteration 19 - deviance = 413.1344 - criterion = 0.0022243
Iteration 20 - deviance = 409.2851 - criterion = 0.00199405
Iteration 21 - deviance = 406.3287 - criterion = 0.001803303
Iteration 22 - deviance = 404.0569 - criterion = 0.001643018
Iteration 23 - deviance = 402.3053 - criterion = 0.001506427
Iteration 24 - deviance = 400.9477 - criterion = 0.001388527
Iteration 25 - deviance = 399.8887 - criterion = 0.001285607
Iteration 26 - deviance = 399.057 - criterion = 0.001194895
Iteration 27 - deviance = 398.3997 - criterion = 0.001114286
Iteration 28 - deviance = 397.8769 - criterion = 0.001042155
Iteration 29 - deviance = 397.459 - criterion = 0.0009772273
Iteration 30 - deviance = 397.1233 - criterion = 0.0009184833
Iteration 31 - deviance = 396.8526 - criterion = 0.0008650971
Iteration 32 - deviance = 396.6336 - criterion = 0.0008163892
Iteration 33 - deviance = 396.456 - criterion = 0.0007717943
Iteration 34 - deviance = 396.3116 - criterion = 0.000730837
Iteration 35 - deviance = 396.1941 - criterion = 0.0006931135
Iteration 36 - deviance = 396.0983 - criterion = 0.0006582784
Iteration 37 - deviance = 396.0201 - criterion = 0.0006260338
Iteration 38 - deviance = 395.9562 - criterion = 0.0005961216
Iteration 39 - deviance = 395.904 - criterion = 0.0005683163
Iteration 40 - deviance = 395.8613 - criterion = 0.0005424203
Iteration 41 - deviance = 395.8264 - criterion = 0.0005182595
Iteration 42 - deviance = 395.7978 - criterion = 0.0004956799
Iteration 43 - deviance = 395.7744 - criterion = 0.0004745448
Iteration 44 - deviance = 395.7553 - criterion = 0.0004547322
Iteration 45 - deviance = 395.7396 - criterion = 0.0004361331
Iteration 46 - deviance = 395.7267 - criterion = 0.0004186497
Iteration 47 - deviance = 395.7162 - criterion = 0.0004021942
Iteration 48 - deviance = 395.7076 - criterion = 0.0003866872
Iteration 49 - deviance = 395.7005 - criterion = 0.0003720569
Iteration 50 - deviance = 395.6948 - criterion = 0.0003582385
Iteration 51 - deviance = 395.69 - criterion = 0.0003451728
Iteration 52 - deviance = 395.6861 - criterion = 0.0003328063
Iteration 53 - deviance = 395.683 - criterion = 0.0003210898
Iteration 54 - deviance = 395.6804 - criterion = 0.0003099788
Iteration 55 - deviance = 395.6782 - criterion = 0.0002994322
Iteration 56 - deviance = 395.6765 - criterion = 0.0002894125
Iteration 57 - deviance = 395.675 - criterion = 0.0002798853
Iteration 58 - deviance = 395.6739 - criterion = 0.0002708189
Iteration 59 - deviance = 395.6729 - criterion = 0.0002621842
Iteration 60 - deviance = 395.6721 - criterion = 0.0002539542
Iteration 61 - deviance = 395.6715 - criterion = 0.0002461042
Iteration 62 - deviance = 395.671 - criterion = 0.0002386112
Iteration 63 - deviance = 395.6705 - criterion = 0.0002314622
Iteration 64 - deviance = 395.6701 - criterion = 0.0002246121
Iteration 65 - deviance = 395.6698 - criterion = 0.0002180685
Iteration 66 - deviance = 395.6696 - criterion = 0.0002118057
Iteration 67 - deviance = 395.6694 - criterion = 0.0002058078
Iteration 68 - deviance = 395.6693 - criterion = 0.0002000602
Iteration 69 - deviance = 395.6691 - criterion = 0.0001945493
Iteration 70 - deviance = 395.669 - criterion = 0.0001892621
Iteration 71 - deviance = 395.6689 - criterion = 0.0001841869
Iteration 72 - deviance = 395.6689 - criterion = 0.0001793124
Iteration 73 - deviance = 395.6688 - criterion = 0.0001746283
Iteration 74 - deviance = 395.6688 - criterion = 0.0001701248
Iteration 75 - deviance = 395.6687 - criterion = 0.0001657927
Iteration 76 - deviance = 395.6687 - criterion = 0.0001616236
Iteration 77 - deviance = 395.6687 - criterion = 0.0001576092
Iteration 78 - deviance = 395.6686 - criterion = 0.0001537422
Iteration 79 - deviance = 395.6686 - criterion = 0.0001500154
Iteration 80 - deviance = 395.6686 - criterion = 0.0001464221
Iteration 81 - deviance = 395.6686 - criterion = 0.000142956
Iteration 82 - deviance = 395.6686 - criterion = 0.0001396112
Iteration 83 - deviance = 395.6686 - criterion = 0.0001363822
Iteration 84 - deviance = 395.6686 - criterion = 0.0001332636
Iteration 85 - deviance = 477.5151 - criterion = 0.001599568
Iteration 86 - deviance = 545.524 - criterion = 0.0001433904
Iteration 87 - deviance = 541.6072 - criterion = 0.0001688429
Iteration 88 - deviance = 531.2352 - criterion = 0.0001711763
Iteration 89 - deviance = 407.89 - criterion = 0.0005527749
Iteration 90 - deviance = 409.2737 - criterion = 0.0001380598
  step size truncated due to possible divergence  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  Inf 
  Stepsize halved - new deviance =  11897.35 
Iteration 91 - deviance = 11897.35 - criterion = 0.7397263
Error in qr.default(x) : NA/NaN/Inf in foreign function call (arg 1)

Traceback() returns:

> traceback()
7: qr.default(x)
6: qr(x)
5: qr(x)
4: matrank(S.k)
3: PQLMQL_innerFit(y.star, X, Z, W, d, groups, offset, method, estimator, 
       control)
2: mmclogit.fitPQLMQL(y = Y, s = s, w = weights, X = XD, Z = ZD, 
       groups = groups, method = method, estimator = estimator, 
       control = control, offset = offset)
1: mblogit(Answer ~ Age_group:Relief:Story, data = long_data2, maxit = 200, 
       method = "PQL", trace = T, random = ~1 | ID)

I then tried changing method to "MQL" but it gets to iteration 99 and returns this error:
Error in qr.default(x) : NA/NaN/Inf in foreign function call (arg 1)

Taking a few of the interactions out does remove this problem of convergence but if possible would like to get an estimate of the full model and all its interactions so I can report the BIC of this full model.

Do you know what I can do here to get the full model to run?

Thanks in advance!
Matthew

mclogit for multinominal logistic regression with fixed effects (conditional multinomial logistic)

Hi, may I ask a question on using this package? I need to run a multinomial logistic regression with fixed effects. It seems that this package is the closest to what I want, but I don't know whether it can do this? if it can, then how to add fixed effects? I have read the package instructions but still don't understand it. Thank you very much!

Error from `predict.mmblogit()`

With mclogit v0.9.4, the following example runs into an error:

data(VA, package = "MASS")
ngrp <- 8L
VA$grp <- gl(n = ngrp,
             k = floor(nrow(VA) / ngrp),
             length = nrow(VA),
             labels = paste0("gr", seq_len(ngrp)))
set.seed(2643)
VA$grp <- sample(VA$grp)

fit_mblogit <- mclogit::mblogit(cell ~ treat + age + Karn,
                                random = ~ 1 | grp,
                                data = VA)

lpreds <- predict(fit_mblogit, newdata = VA)

The last line throws

Error in `[.data.frame`(mf, group.labels) : undefined columns selected

which I guess is due to the changes in predict.mmblogit() from commit 371692d.

Incorrect predictions from mblogit when predictors are scaled

I have found that when an mblogit model contains scaled predictor(s), then subsequent predictions are incorrect. Here is an example comparing predictions from equivalent models using nnet::multinom and mclogit::mblogit:

library(MASS)
library(nnet)
library(mclogit)
#> Loading required package: Matrix

# silly modification to housing data
housing = transform(MASS::housing, x = 13 + as.numeric(Infl))

# x values for testing with prediction
newx = data.frame(x = 10:12)

# Model with multinom:
h.mult = multinom(Sat ~ scale(x), weights = Freq, data = housing)
#> # weights:  9 (4 variable)
#> initial  value 1846.767257 
#> final  value 1772.112751 
#> converged
predict(h.mult, newdata = newx, type = "probs")
#>         Low    Medium       High
#> 1 0.8326790 0.1448212 0.02249979
#> 2 0.7703353 0.1843759 0.04528885
#> 3 0.6862040 0.2260201 0.08777590

# Model with mblogit
h.mbl = mblogit(Sat ~ scale(x), weights = Freq, data = housing)
#> 
#> Iteration 1 - Deviance = 3557.833
#> Iteration 2 - Deviance = 3544.227
#> Iteration 3 - Deviance = 3544.226
#> Iteration 4 - Deviance = 3544.226
#> converged
predict(h.mbl, newdata = newx, type = "response")
#>            Low    Medium      High
#> [1,] 0.4245383 0.2802862 0.2951755
#> [2,] 0.3148686 0.2702901 0.4148413
#> [3,] 0.2167930 0.2419704 0.5412366

Note that the predictions differ substantially. This is in part because the scaling information is not available. Normally, this information is in the predvars attribute of the terms component:

attr(h.mult$terms, "predvars")
#> list(Sat, scale(x, center = 15, scale = 0.822226451793038))
attr(h.mbl$terms, "predvars")
#> NULL

I tried tricking it into doing the right thing by copying this attribute from the other model. However, this does not change the predictions:

attr(h.mbl$terms, "predvars") = attr(h.mult$terms, "predvars")
predict(h.mbl, newdata = newx, type = "response")
#>            Low    Medium      High
#> [1,] 0.4245383 0.2802862 0.2951755
#> [2,] 0.3148686 0.2702901 0.4148413
#> [3,] 0.2167930 0.2419704 0.5412366

^{Created on 2021-02-16 by the reprex package (v0.3.0)}

A few notes:

Presumably, similar errors will occur with other predictor functions, e.g., poly(), where we need the transformation parameters based on the original data.
Two changes are needed to set this right:
a) Make sure the $terms component of the returned model is complete. Normally, it can be obtained as an attribute of the model matrix
b) The predict.mblogit method seems to build its model matrix X from the model formula (In see rhs <- object$formula[-2] in the code), rather than the $terms component of the object.

Problem in a model with two nested variables

Hello Prof. Martin Elff,

I'm trying to fit a multinomial logit model with two random parameters, since my data set has individuals that are enrolled in specific institutions and programs within institutions. The mblogit function allow to do that? If so, how can I specify the formula? I tried to follow the logic from the lme4 package, but I get the following results:

# A tibble: 10 x 7
         id_ institution sex   age_group      secondary_school dropout  program_id
       <dbl> <fct>       <fct> <fct>          <fct>            <fct>    <fct>     
 1 1.98e-312 3849        M     19 to 24 years Public school    Dropout  384918424 
 2 1.98e-312 17          M     19 to 24 years Public school    Enrolled 171439    
 3 1.98e-312 600         M     Above 30 years Public school    Enrolled 6001103914
 4 1.98e-312 694         F     Up to 18       Private school   Enrolled 69415839  
 5 1.98e-312 573         F     19 to 24 years Public school    Stopout 573116884 
 6 1.98e-312 694         M     Above 30 years Public school    Enrolled 69452141  
 7 1.98e-312 5           M     Up to 18       Private school   Enrolled 5300518   
 8 1.98e-312 601         M     19 to 24 years Public school    Dropout  6011327404
 9 1.98e-312 2564        F     19 to 24 years Private school   Enrolled 2564118910
10 1.98e-312 578         F     19 to 24 years Public school    Enrolled 57813288

mblogit(formula = dropout ~ sex + age_group + secondary_school,
random = ~ 1 | institution/program_id,
data = stud.db)

Iteration 1 - Deviance = 79497.12
Iteration 2 - Deviance = 75292.89
Iteration 3 - Deviance = 74943.24
Iteration 4 - Deviance = 74931.05
Iteration 5 - Deviance = 74931
Iteration 6 - Deviance = 74931
converged

Error in Z[[k]] : subscript out of bounds

The class_id variables is a index that combines the id from the institution and the id from the program within the institution (because some programs appears in more than one institution).

Thanks a lot.

Best,

`tmp`[[k]] : subscript out of bounds

Dear Dr. Elff,

I was trying to run function mblogit today (version 0.8.7.1) for a multinomial mixed logistic regression model. Here is the code:

mblogit(value ~ Scan_type , random = ~ 1|(Rater2 + Subject_ID), data = final_df)

This model works fine with 1 random intercept (either Rater2 or Subject_ID), but when consider both of them, I have this error after 6 iterations (models converge too!): *tmp*[[k]] : subscript out of bounds.

Could you please let me know what I could do to get around this problem?

Thanks,
Quy

Error when using nested random terms

Hello and thank you for developing this package. I am running into a persistent error message while trying to construct a mixed effect, baseline-category multinomial logit and hope you can provide some insight.

I am trying to construct a model which describes how different factors influence the likelihood of several behavioral states among two species of shark. 19 individual animals, comprising two distinct species, were each tracked for ~100 days each and a different behavioral state was identified for each day data were collected. I would like to code individual shark as a categorical random effect variable (with 19 levels) within species, a categorical fixed effect variable (with 2 levels).

With this general idea, the code that I am currently trying to run is:
mclogit::mblogit(cluster ~ 1, random = ~1|species/individual, data = df)

Which produces the error message:

Error in solve.default(X[[i]], ...) : 'a' (6 x 1) must be square

This error message is only produced in models where 'species' is included as a random effect and not when species is included as a fixed effect or individual is included along. This suggests that the error is in how the function is processing the 'species' variable specifically.

Here is a dummy sample of the data on which I am trying to train the model:

df <- data.frame(Date = c("2015-11-25", "2016-01-24", "2016-02-27", "2016-03-27", "2017-12-02", "2017-12-06", "2015-10-30", "2015-10-31", "2009-01-28", "2009-02-04", "2009-02-15", "2009-02-20", "2009-02-21", "2009-02-23"), cluster = factor(c(3,3,4,6,3,1,3,2, 2, 4, 3, 3, 3, 4)),species = factor(c("I.oxyrinchus", "I.oxyrinchus", "I.oxyrinchus", "I.oxyrinchus", "P.glauca", "P.glauca", "P.glauca", "P.glauca", "I.oxyrinchus", "I.oxyrinchus", "I.oxyrinchus", "P.glauca", "P.glauca", "P.glauca")),individual = factor(c("141257", "141257", "141254", "141254", "141256", "141256", "141255", "141255", "78680", "78682", "78683", "78680", "78682", "78683")))

Similar to the whole dataset, this model runs this code:
mclogit::mblogit(cluster ~ 1, random = ~1|individual, data = df)
But fails to run these:
mclogit::mblogit(cluster ~ 1, random = ~1|species/individual, data = df)
mclogit::mblogit(cluster ~ 1, random = ~1|species, data = df)

Interestingly, with the factor variable 'individual' has fewer distinct levels, it produces the same error message. Perhaps the breakpoint exists in how the function handles random variables with few discrete levels.

After following the error traceback the breakpoint seems to originate from:

lapply(Phi.start, solve) > PQLMQL_innerFit(y.star, X, Z, W, d, offset, method, estimator, control) > mmclogit.fitPQLMQL(y = Y, s = s, w = weights, X = XD, Z = ZD, d = d, method = method, estimator = estimator, control = control, offset = offset)

Based on a little digging, it seems that the error may lie in the object Z.

Can you please provide any more information about what is producing this error message and how I might work around it?
Thank you very much for your time and assistance.

Cheers

Error `'a' (<dimension1> x 1) must be square`

For some models fitted by mclogit::mblogit(), I get the error

Error in solve.default(X[[i]], ...) : 'a' (<dimension1> x 1) must be square

thrown by mclogit:::PQLMQL_innerFit(). (The <dimension1> can be 3, for example.)

I think the problematic line is

mclogit/pkg/R/mmclogit-fitPQLMQL.R

Line 271 in 7f47d53

S.k <- diag(x=S.k,nrow=d)

because, if S.k is a matrix there, then diag() returns a vector containing the diagonal elements. Perhaps that line was supposed to read

S.k <- diag(x=diag(S.k),nrow=d)

instead? But I don't know anything about the mathematical details of this line, so my guess might be completely wrong.

And I'm not sure, but perhaps this is the same bug as in #21?

Note that this issue occurs with both, the current CRAN version 0.8.7.3 as well as with the current GitHub version at commit 7f47d53 installed via devtools::install_github("melff/mclogit",subdir="pkg").

non-convergence and errors with random effects models using mblogit

Hi. I was using the CRAN version of package mclogit and none of my random effects models using the mblogit command converged, whether or not I used the control argument. I read an online post from another user who encountered this and said it was a package error that had been fixed in the version currently on Github, so I installed that version. The random effects model still fails to converge if I do not use the control argument and the summary warns, "In sqrt(diag(vcov.phi)) : NaNs produced." I understand my data may simply be unable to support a random effects model. However, I tried, but failed, to implement the control argument in the Github version of the package.

Here is the command and error message when I use control = mmclogit.control

MBLrW <- mblogit(aBehavior ~ Trt15_25_5min, random=~1|Deployment, data=HO,
control=mmclogit.control(epsilon = 1e-12, maxit = 100, trace=TRUE,
trace.inner=FALSE, avoid.increase = FALSE,break.on.increase = FALSE,
break.on.infinite = FALSE, break.on.negative = FALSE))

Error in mmclogit.control(epsilon = 1e-12, maxit = 100, trace = TRUE, :

could not find function "mmclogit.control"

Here is the command and error message when I use control = mclogit.control

MBLrW <- mblogit(aBehavior ~ Trt15_25_5min, random=~1|Deployment, data=HO,
control=mclogit.control(epsilon = 1e-12, maxit = 100, trace=TRUE))

Error in mmclogit.fitPQLMQL(y = Y, s = s, w = weights, X = XD, Z = ZD, : object 'step.truncated' not found

In addition: Warning messages:

1: In mmclogit.fitPQLMQL(y = Y, s = s, w = weights, X = XD, Z = ZD, : Numeric problems in inner iteration, bailing out

2: Algorithm did not converge

Anything you could do to help me would be wonderful. Thank you.

Problem with random effect: In sqrt(diag(vcov.phi)) : NaNs produced

Hello Professor,

I'm using your package for a meta-analysis.
And I don't understand why I get this error message with my random effect "Ref": Warning message: In sqrt(diag(vcov.phi)) : NaNs produced.
Is there a solution to fix this? Or is it just due to the heterogeneity of my dataset (with some "Ref" where I have 100 observations and others only 1)?
How do I interpret my model when it displays this error?

Call:
mblogit(formula = y ~ x1 + x2 + scale(x3)
    data = data_trends, random = ~1 | Ref, method = "PQL", estimator = "ML", 
    control = mmclogit.control(epsilon = 1e-08, maxit = 1000, 
        trace = TRUE, trace.inner = FALSE, avoid.increase = FALSE, 
        break.on.increase = FALSE, break.on.infinite = T, break.on.negative = FALSE))

Equation for B vs A
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)           1.26436    0.21421   5.902 3.58e-09 ***
x1Def                -1.35557    0.14634  -9.263  < 2e-16 ***
x2Org                -0.08302    0.27047  -0.307  0.75888    
scale(x3)            -0.41593    0.12978  -3.205  0.00135 ** 

Equation for C vs A
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)            1.8383     0.2028   9.064  < 2e-16 ***
x1Def                 -1.3911     0.1370 -10.152  < 2e-16 ***
x2Org                 -0.3747     0.2565  -1.461    0.144    
scale(x3)             -0.5783     0.1202  -4.809 1.52e-06 ***

Equation for D vs A
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)           -1.2673     0.2704  -4.686 2.78e-06 ***
x1Def                 -0.3602     0.2280  -1.580   0.1142    
x2Org                  0.5224     0.2883   1.812   0.0700 .  
scale(x3)             -0.2753     0.1365  -2.017   0.0437 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Co-)Variances:

Grouping level: Ref 
           Estimate                                    Std.Err.                        
      B~1        C~1        D~1           B~1        C~1            D~1    
B~1   1.7604                         B~1  0.09226                        
C~1   0.7066     1.6297              C~1  0.25405    0.17952             
D~1   0.2596     0.4086     0.9585   D~1  0.11592        NaN        NaN  

Approximate residual deviance: 5368 
Number of Fisher scoring iterations:  7
Number of observations
  Groups by Ref: 156
  Individual observations:  2595

Thanks in advance

Aggregate logit

Although the mclogit function can calculate logit models for not only count (the Molel 1 below) but also share/ratio (the Model 2 below), the detail for the model 2 is currently not described in the manual. So, I would like to know whether the model 2 is correspond to the aggregate logit model or not.

library(mclogit)
data(Transport)
### Model 1: for count
mclogit(cbind(resp, suburb)~distance+cost,data=Transport) 
### Model 2: for share/ratio 
mclogit(cbind(prop.true, suburb)~distance+cost,data=Transport)

Best wishes and thank you so much for the very flexible package!

Feature request: allow some form of regularisation in multinomial fits?

More as a feature request: do you think it might be possible to allow some form of regularisation in mblogit multinomial fits?
When one is dealing with sparse data, one now often encounters problems of complete separation, resulting in very large errors. This could be resolved if some form of regularisation could be added. I understand that for logistic regression, a ridge (L2 norm) penalty on the coefficients can be added just by changing the IRLS algorithm (https://bwlewis.github.io/GLM/) so that in the least square regression step of the IRLS algorithm it uses a ridge regression instead of a regular OLS regression. This can be done just by augmenting the covariate matrix X with a p x p covariate matrix with p = nr of columns of the covariate matrix and sqrt(lambdas) along the diagonal and augmenting the outcome variable with p zeros (with lambdas=a vector of L2 norm penalty coefficients=lambda*penalty weights, with vector penalty weights typically being set to 0 for the intercept and to 1 for the rest, or in case of adaptive ridge regression to 1/(abs(coef)) or 1/(coef^2), where coef are the coefficients obtained e.g. in an earlier ridge regression) (see https://stats.stackexchange.com/questions/69205/how-to-derive-the-ridge-regression-solution and https://stats.stackexchange.com/questions/137057/ridge-penalized-glms-using-row-augmentation). I presume the same recipe could be applied for a multinomial regression, in which case incorporating such a feature might not be too hard - in that case it could just be a matter of changing line 37 in https://github.com/melff/mclogit/blob/master/pkg/R/mclogit-fit.R to a weighted least square ridge regression.

For example, in a multinomial model looking at SARS-CoV2 lineage dynamics through time in different countries & continents
lineage ~ date+date:continent+date:country+country
it would be nice if I could penalize just the date:country logistic slopes, and leave all the intercepts & the per-continent logistic slopes unpenalized, in order to shrink the per-country logistic slopes a little towards the ones of the continent where the country is in.

Do you think such a feature could be incorporated in the future?

I saw that package MGLM tried to incorporate LASSO penalties for multinomial models in their function
https://www.rdocumentation.org/packages/MGLM/versions/0.2.1/topics/MGLMsparsereg
but I had no luck getting it to work (and it also didn't provide a vcov method).
And glmnet of course allows for elastic net regularised multinomial models (but that one doesn't return the variance-covariance matrix and doesn't have a formula interface).

In terms of regularised (L1/LASSO or L2/ridge penalized) multinomial models I also saw this one,
https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1424458. Though I believe the variance-covariance matrix only has a closed-form solution for the ridge penalized case.

Error with "predict" command in newest versions

Hello Prof. Elff,
I'm using the mclogit command for an analysis and encountered an error with the post estimation command "predict" when using a more recent version of the package. The error reads:

Error in .TM.repl.i.mat(as(x, "TsparseMatrix"), j=j, value=value) : column indices must be <= ncol(.) which is 1417.

I'm able to get the predicted choice probabilities with version 0.8.5.1, but the error with the newer version makes me think that there may be a bug. Thanks for all of your work on this package; it's been incredibly useful for me.
Cheers,
Kaylyn

mblogit model comparison with anova does not calculate degrees of freedom from random parameters

Thanks for the great package!

I ran into a minor problem when comparing models. I tried to run two models, one with random intercepts and other with random intercepts and slopes. Model comparison with "anova" did calculate the deviance difference correctly, but there seems to be something wrong with degrees of freedom. It seems to only include the fixed formula in this calculation. Below is an example

library(mclogit)

set.seed(21385)
n=2000

d<-data.frame(
  DV=as.factor(sample(c("A","B","C"),n,replace=T)),
  gender=sample(c(-0.5,0.5),n,replace=T),
  cntry=sample(letters[1:10],n,replace=T))



#with random intercept
fit.1<-mblogit(DV~gender,random=~1|cntry,data=d)
summary(fit.1)

#with random slope
fit.2<-mblogit(DV~gender,random=~1+gender|cntry,data=d)
summary(fit.2)


#anova reads same number of degrees of freedom for both models
anova(fit.2,fit.1,test="Chisq")

Incorrect result by `predict.mmblogit()` in case of multiple `random` formulas

I'm not 100% sure, but I think line

mclogit/pkg/R/mblogit.R

Line 872 in dadb418

ZD <- blockMatrix(ZD)

needs to be removed. Consider the following reprex:

data(VA, package = "MASS")
ngrp <- 8L
VA$grp <- gl(n = ngrp,
             k = floor(nrow(VA) / ngrp),
             length = nrow(VA),
             labels = paste0("gr", seq_len(ngrp)))
set.seed(2643)
VA$grp <- sample(VA$grp)

nsbj <- 13L
VA$sbj <- gl(n = nsbj,
             k = floor(nrow(VA) / nsbj),
             length = nrow(VA),
             labels = paste0("subj", seq_len(nsbj)))
VA$sbj <- sample(VA$sbj)

fit_multi <- mclogit::mblogit(cell ~ treat + age + Karn,
                              random = list(~ 1 | sbj,
                                            ~ 1 | grp),
                              data = VA)

VA_new <- head(VA, 1)
VA_new$sbj <- "subj1"
VA_new$grp <- "gr1"
lpreds_multi <- predict(fit_multi, newdata = VA_new)

# Manual check:
VA_new$treat2 <- as.numeric(VA_new$treat == "2")
fixnms_b <- c("treat2", "age", "Karn")
nlats <- nlevels(VA$cell) - 1L
mfixef <- fit_multi$coefmat
mranef <- setNames(fit_multi$random.effects, names(fit_multi$groups))
# Actual result:
lpreds_multi_man <- matrix(mfixef[, "(Intercept)"],
                           nrow = nrow(VA_new),
                           ncol = nlats, byrow = TRUE) +
  as.matrix(VA_new[, fixnms_b, drop = FALSE]) %*% t(mfixef[, fixnms_b]) +
  t(mranef$sbj[1:3, ])
all.equal(unname(lpreds_multi), unname(lpreds_multi_man), tolerance = 1e-15)
## --> Gives TRUE.
# Expected result:
lpreds_multi_man_expect <- lpreds_multi_man +
  t(mranef$grp[1:3, ])
all.equal(unname(lpreds_multi), unname(lpreds_multi_man_expect), tolerance = 1e-15)
## --> Gives "Mean relative difference: 0.2396385".

# After commenting line `ZD <- blockMatrix(ZD)` in mclogit:::predict.mmblogit() and
# running
lpreds_multi <- predict(fit_multi, newdata = VA_new)
# I get
all.equal(unname(lpreds_multi), unname(lpreds_multi_man), tolerance = 1e-15)
## --> Gives "Mean relative difference: 0.2759492".
# and
all.equal(unname(lpreds_multi), unname(lpreds_multi_man_expect), tolerance = 1e-15)
## --> Gives TRUE.
# as expected.

EDIT: I was using mclogit version 0.9.4.1 (at commit dadb418) for this.

Possible typo in mblogit code

In the code in mblogit for dealing with the groups argument, there seems to be a typo on line 124:

gf <- as.formula(rf)

As far as I can tell, the variable rf should never be defined at this point in the code (since it is only defined in a mutually exclusive branch of the if/else), so this line will always throw an error. Presumably the correct line is:

gf <- as.formula(gf)