melff / mclogit Goto Github PK
View Code? Open in Web Editor NEWmclogit: Multinomial Logit Models, with or without Random Effects or Overdispersion
Home Page: http://melff.github.io/mclogit/
mclogit: Multinomial Logit Models, with or without Random Effects or Overdispersion
Home Page: http://melff.github.io/mclogit/
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
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.
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?
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))
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))
Anything you could do to help me would be wonderful. Thank you.
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.
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")
(om = mblogit(y~x1+x2, random=~1|country, data=fake))
(om = mblogit(y~x1+x2, random=~1|county/country, data=fake))
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.
I'm not 100% sure, but I think line
Line 872 in dadb418
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.
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:
poly()
, where we need the transformation parameters based on the original data.$terms
component of the returned model is complete. Normally, it can be obtained as an attribute of the model matrixpredict.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.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,
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.
Calling summary
on an mblogit
model fit does odd things that depend on the level names of the response. For example:
set.seed(123)
d = data.frame(R=factor(sample(c("A","B","C"),30,TRUE)), X1=rnorm(30), C1 = factor(sample(letters[1:2],30,TRUE)))
m1 = mblogit(R~X1+C1, data=d, trace=FALSE)
levels(d$R) = c("A","BC","C")
m2 = mblogit(R~X1+C1, data=d, trace=FALSE)
m1
and m2
look correct when printed:
> m1
Call: mblogit(formula = R ~ X1 + C1, data = d, trace = FALSE)
Coefficients:
Predictors
Response categories (Intercept) X1 C1b
B/A 1.2430 0.5280 -1.6924
C/A 0.7349 0.4068 -0.4894
Null Deviance: 65.92
Residual Deviance: 61.44
> m2
Call: mblogit(formula = R ~ X1 + C1, data = d, trace = FALSE)
Coefficients:
Predictors
Response categories (Intercept) X1 C1b
BC/A 1.2430 0.5280 -1.6924
C/A 0.7349 0.4068 -0.4894
Null Deviance: 65.92
Residual Deviance: 61.44
but summary
prints additional terms
> summary(m1)
[cut]
Equation for B vs A:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.2430 0.7560 1.644 0.100
X1 0.5280 0.5496 0.961 0.337
C1b -1.6924 1.0415 -1.625 0.104
Equation for C vs A:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.7349 0.8021 0.916 0.360
X1 0.4068 0.5293 0.769 0.442
C1b -0.4894 0.9840 -0.497 0.619
[cut]
> summary(m2)
[cut]
Equation for BC vs A:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.2430 0.7560 1.644 0.100
X1 0.5280 0.5496 0.961 0.337
C1b -1.6924 1.0415 -1.625 0.104
Equation for C vs A:
Estimate Std. Error z value Pr(>|z|)
B(Intercept) 1.2430 0.7560 1.644 0.100
(Intercept) 0.7349 0.8021 0.916 0.360
BX1 0.5280 0.5496 0.961 0.337
X1 0.4068 0.5293 0.769 0.442
BC1b -1.6924 1.0415 -1.625 0.104
C1b -0.4894 0.9840 -0.497 0.619
[cut]
These additional terms (B(Intercept)
, BX1
and BC1b
) appear to come from the Equation for BC vs A, and this happens when a factor level has another factor level as part of its name - in this case "BC" ends with "C". Originally this happened when I was modelling with factors including "Medicine" and "Emergency Medicine", and I tracked it down to the minimal example above.
My workaround is currently to change the response level names so this doesn't happen. It doesn't seem to affect the actual model fits, just the output from summary.
This is using 0.9.6 version of mclogit, from CRAN.
Hi,
I've been trying to set up a multinomial model with random effects using mclogit, but I get stuck because no matter which dataset or formula I use, I always get the same error in R "Error in formula[[2]][[1]] : object of type 'symbol' is not subsettable". I've looked for a solution around but could not find one, any ideas how to fix this?
Any help would be much appreciated! Thanks!
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
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")
.
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
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.
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!
I'm trying to fit a multinomial response mixed effects model on a largish (~150,000 row) data set and this package seems to be the only that can do it in reasonable time. But now I'm struggling to compare the model outputs with outputs from other mixed effects modelling packages, specifically lme4.
With an lme4 fitted model I can do ranef(m)
and get the random effects, with the conditional variances matrix too. This can be fed into dotplot
to get a plot of "confidence intervals" for the random effects. See help(lme4::ranef)
for an example.
What I want to do is make these plots from mblogit
output. It seems the $random.effects
attribute has the random effects values, but I'm struggling to find or compute the equivalent conditional variances matrix.
My current working assumption is that the following two binomial models have the same model specification:
library(lme4)
library(HSAUR3) # toenail data
library(mclogit)
gl <- glmer(outcome~visit+(1|treatment), data=toenail, family=binomial)
mb = mblogit(outcome~visit, random=~1|treatment, data=toenail)
and if I can replicate ranef(gl)
from mb
then I'll be most of the way there. There are differences in these fitted models but I'm not sure if that's down to a completely different model specification or just a different optimiser stopping in essentially the same flat point of likelihood space.
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)
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)
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
Hello- I have two loosely related question maybe you can help. A bit of background:
I have data about the forms that a parasite can take in infected mice, so each datapoint is a parasite. Parasites are collected within mice, and mice within experiments. Each mouse is infected with one of four lines of parasites. The question of interest is the effect of parasite line on the proportion of the three forms accounting for the effect of mouse and experiments. In particular, I'm interested in the difference between the WT line and each of the three KO lines. This is a summary of number of mice, lines and experiments:
Line
Experiment KO_1 KO_3 KO_7 WT
A 0 3 0 3
B 0 3 0 2
C 0 0 3 3
D 3 0 3 3
E 3 0 3 3
The model I have is:
bfitu <- mblogit(Form ~ Line, random=list(~1|Experiment/mouse_id), data=lldat)
emmeans::emmeans(bfitu, ~ Line | Form)
This gives sensible results but I get the warnings below. Should I be worried about them?
Warning: Inner iterations did not coverge - nlminb message: false convergence (8)
Iteration 1 - deviance = 12516.11 - criterion = 0.9603497Warning: Inner iterations did not coverge - nlminb message: false convergence (8)
Iteration 2 - deviance = 12430.75 - criterion = 0.05545146Warning: Inner iterations did not coverge - nlminb message: false convergence (8)
Iteration 3 - deviance = 12421.89 - criterion = 0.01194718Warning: Inner iterations did not coverge - nlminb message: false convergence (8)
Iteration 4 - deviance = 12420.05 - criterion = 0.001537364Warning: Inner iterations did not coverge - nlminb message: false convergence (8)
Iteration 5 - deviance = 12419.78 - criterion = 2.682554e-05Warning: Inner iterations did not coverge - nlminb message: false convergence (8)
Iteration 6 - deviance = 12419.75 - criterion = 6.59256e-09
converged
Second question: I'm not entirely sure that a more appropriate model would be:
bfitu2 <- mblogit(Form ~ Line, random=list(~1|Experiment, ~1|mouse_id), data=lldat)
any advise on this would be much appreciated.
The dataset is attached (data.tsv.gz)
Thanks for the great package!
sessionInfo()
R version 4.3.2 (2023-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS
Matrix products: default
BLAS/LAPACK: /home/dario/miniforge3/envs/20240220_multinomial/lib/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
[5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
[7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
time zone: Europe/London
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggbeeswarm_0.7.2 emmeans_1.10.0 mclogit_0.9.6 Matrix_1.6-5
[5] nnet_7.3-19 ggplot2_3.4.4 data.table_1.14.10
loaded via a namespace (and not attached):
[1] gtable_0.3.4 jsonlite_1.8.8 dplyr_1.1.4 compiler_4.3.2
[5] tidyselect_1.2.0 scales_1.3.0 yaml_2.3.8 lattice_0.22-5
[9] coda_0.19-4.1 R6_2.5.1 generics_0.1.3 MASS_7.3-60
[13] tibble_3.2.1 munsell_0.5.0 pillar_1.9.0 rlang_1.1.3
[17] utf8_1.2.4 estimability_1.4.1 cli_3.6.2 withr_3.0.0
[21] magrittr_2.0.3 grid_4.3.2 mvtnorm_1.2-4 xtable_1.8-4
[25] memisc_0.99.31.7 beeswarm_0.4.0 lifecycle_1.0.4 vipor_0.4.7
[29] vctrs_0.6.5 glue_1.7.0 farver_2.1.1 fansi_1.0.6
[33] colorspace_2.1-0 tools_4.3.2 pkgconfig_2.0.3
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!
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!
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")
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
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.
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
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
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.