I am running multinomial regression using nnet::multinom
for a dataset where response and predictors are all factors. A sample of my data:
df2 = structure(list(SOURC = c("dsint", "dsint", "circm", "dsint",
"dsint", "circm", "dsint", "circm", "circm", "dsint", "dsint",
"circm", "circm", "dsint", "circm", "circm", "circm", "circm",
"dsint", "rul", "circm", "dsint", "dsint", "dsint", "dsint",
"dsint", "dsint", "circm", "dsint", "circm", "dsint", "dsint",
"circm", "circm", "circm", "cond", "circm", "dsint", "dsint",
"circm", "circm", "dsint", "dsint", "circm", "dsint", "dsint",
"circm", "circm", "circm", "circm", "dsint", "dsint", "circm",
"rul", "circm", "dsint", "circm", "dsint", "circm", "circm",
"dsint", "circm", "circm", "dsint", "circm", "dsint", "dsint",
"circm", "dsint", "dsint", "dsint", "dsint", "rul", "cond", "circm",
"cond", "circm", "dsint", "dsint", "dsint", "circm", "circm",
"dsint", "dsint", "dsint", "dsint", "dsint", "dsint", "circm",
"circm", "dsint", "dsint", "dsint", "dsint", "dsint", "dsint",
"circm", "circm", "dsint", "dsint"), SITN = c("dynamic", "dynamic",
"dynamic", "stative", "dynamic", "dynamic", "dynamic", "dynamic",
"dynamic", "dynamic", "dynamic", "dynamic", "dynamic", "dynamic",
"dynamic", "dynamic", "dynamic", "dynamic", "dynamic", "dynamic",
"dynamic", "stative", "dynamic", "stative", "dynamic", "dynamic",
"dynamic", "dynamic", "dynamic", "dynamic", "stative", "stative",
"dynamic", "dynamic", "dynamic", "stative", "dynamic", "dynamic",
"dynamic", "dynamic", "dynamic", "dynamic", "dynamic", "stative",
"stative", "stative", "dynamic", "dynamic", "dynamic", "dynamic",
"dynamic", "stative", "stative", "dynamic", "stative", "dynamic",
"dynamic", "dynamic", "dynamic", "dynamic", "dynamic", "dynamic",
"dynamic", "stative", "dynamic", "dynamic", "dynamic", "stative",
"dynamic", "stative", "dynamic", "dynamic", "dynamic", "dynamic",
"stative", "dynamic", "dynamic", "stative", "dynamic", "stative",
"dynamic", "dynamic", "dynamic", "dynamic", "dynamic", "stative",
"dynamic", "dynamic", "dynamic", "dynamic", "dynamic", "stative",
"stative", "stative", "stative", "dynamic", "dynamic", "dynamic",
"stative", "dynamic"), GENR = c("specific", "specific", "generic",
"specific", "specific", "specific", "generic", "generic", "specific",
"specific", "specific", "specific", "specific", "specific", "generic",
"specific", "generic", "specific", "specific", "generic", "generic",
"specific", "specific", "specific", "specific", "specific", "generic",
"generic", "specific", "generic", "generic", "specific", "specific",
"specific", "specific", "specific", "specific", "specific", "generic",
"generic", "specific", "specific", "specific", "specific", "specific",
"generic", "specific", "specific", "specific", "generic", "specific",
"generic", "generic", "generic", "specific", "specific", "specific",
"specific", "generic", "specific", "specific", "specific", "specific",
"generic", "specific", "generic", "generic", "specific", "specific",
"specific", "specific", "specific", "generic", "generic", "specific",
"specific", "generic", "generic", "specific", "specific", "specific",
"specific", "specific", "specific", "specific", "generic", "specific",
"generic", "specific", "specific", "generic", "generic", "specific",
"specific", "specific", "generic", "generic", "specific", "specific",
"specific"), PRSN = c("1st", "3rd", "2nd", "1st", "1st", "3rd",
"2nd", "2nd", "3rd", "3rd", "3rd", "1st", "3rd", "3rd", "1st",
"3rd", "3rd", "1st", "3rd", "3rd", "2nd", "3rd", "3rd", "2nd",
"3rd", "3rd", "3rd", "2nd", "1st", "2nd", "2nd", "3rd", "3rd",
"1st", "1st", "3rd", "1st", "3rd", "2nd", "2nd", "1st", "1st",
"3rd", "1st", "3rd", "3rd", "2nd", "1st", "1st", "1st", "1st",
"3rd", "3rd", "3rd", "3rd", "1st", "3rd", "3rd", "3rd", "1st",
"1st", "2nd", "1st", "1st", "3rd", "1st", "2nd", "3rd", "3rd",
"1st", "3rd", "3rd", "2nd", "3rd", "1st", "2nd", "2nd", "3rd",
"3rd", "3rd", "1st", "2nd", "2nd", "1st", "3rd", "3rd", "3rd",
"3rd", "1st", "3rd", "3rd", "3rd", "3rd", "1st", "1st", "2nd",
"2nd", "3rd", "1st", "2nd"), ANIM = c("animate", "inanimate",
"animate", "animate", "animate", "animate", "animate", "animate",
"inanimate", "animate", "inanimate", "animate", "inanimate",
"animate", "animate", "animate", "animate", "animate", "inanimate",
"inanimate", "animate", "animate", "animate", "animate", "animate",
"animate", "inanimate", "animate", "animate", "animate", "animate",
"inanimate", "animate", "animate", "animate", "inanimate", "animate",
"animate", "animate", "animate", "animate", "animate", "animate",
"animate", "animate", "animate", "animate", "animate", "animate",
"animate", "animate", "inanimate", "inanimate", "animate", "animate",
"animate", "inanimate", "animate", "animate", "animate", "animate",
"animate", "animate", "animate", "inanimate", "animate", "animate",
"animate", "inanimate", "animate", "animate", "animate", "animate",
"animate", "animate", "animate", "animate", "inanimate", "inanimate",
"animate", "animate", "animate", "animate", "animate", "inanimate",
"animate", "animate", "animate", "animate", "inanimate", "animate",
"animate", "inanimate", "animate", "animate", "animate", "animate",
"animate", "animate", "animate"), AGENT = c("yes", "no", "yes",
"yes", "yes", "yes", "yes", "yes", "no", "yes", "no", "yes",
"no", "yes", "yes", "yes", "yes", "yes", "no", "no", "yes", "no",
"yes", "yes", "yes", "yes", "no", "yes", "yes", "yes", "no",
"no", "yes", "yes", "yes", "no", "yes", "yes", "yes", "yes",
"yes", "yes", "no", "no", "no", "no", "yes", "yes", "yes", "yes",
"yes", "no", "no", "yes", "no", "yes", "no", "yes", "yes", "yes",
"yes", "yes", "yes", "yes", "no", "yes", "yes", "no", "no", "no",
"no", "yes", "yes", "yes", "no", "yes", "yes", "no", "no", "no",
"yes", "yes", "yes", "no", "no", "yes", "yes", "yes", "yes",
"no", "yes", "no", "no", "no", "yes", "yes", "yes", "yes", "no",
"yes"), CTGRY = c("cmt", "cmt", "wbf", "twt", "wbf", "wbf", "wbs",
"wbf", "wbf", "cmt", "cmt", "twt", "wbf", "cmt", "wbs", "cmt",
"cmt", "twt", "wbf", "wbs", "wbs", "cmt", "cmt", "twt", "cmt",
"cmt", "cmt", "wbs", "cmt", "wbs", "cmt", "wbf", "wbf", "wbf",
"wbs", "cmt", "twt", "cmt", "wbs", "wbs", "wbs", "cmt", "cmt",
"wbf", "wbf", "cmt", "wbf", "cmt", "cmt", "wbs", "twt", "cmt",
"wbf", "wbf", "wbf", "wbf", "wbf", "cmt", "cmt", "twt", "cmt",
"wbf", "cmt", "cmt", "wbf", "cmt", "wbs", "cmt", "wbs", "twt",
"cmt", "cmt", "wbs", "wbs", "twt", "wbf", "wbs", "wbf", "cmt",
"wbf", "wbf", "twt", "twt", "cmt", "cmt", "cmt", "wbf", "cmt",
"twt", "wbs", "wbf", "wbf", "cmt", "cmt", "cmt", "twt", "wbs",
"cmt", "wbf", "twt"), CNTRY = c("in", "in", "us", "in", "us",
"in", "in", "us", "us", "us", "in", "in", "in", "us", "in", "in",
"in", "in", "in", "in", "in", "in", "in", "us", "us", "in", "in",
"in", "in", "us", "us", "in", "us", "us", "in", "us", "us", "in",
"in", "in", "in", "in", "us", "in", "in", "us", "in", "in", "in",
"us", "us", "in", "us", "in", "in", "in", "us", "us", "in", "us",
"in", "in", "us", "us", "in", "us", "in", "in", "in", "in", "us",
"in", "us", "us", "in", "in", "in", "in", "in", "in", "in", "us",
"in", "us", "in", "in", "us", "in", "in", "in", "us", "in", "in",
"in", "us", "us", "in", "in", "us", "us"), VERB = c("ndt", "hvt",
"hvt", "ndt", "hvt", "hvt", "ndt", "ndt", "ndt", "mst", "ndt",
"ndt", "ndt", "ndt", "mst", "hvt", "hvt", "ndt", "hvt", "mst",
"mst", "mst", "mst", "hvt", "hvt", "mst", "ndt", "ndt", "hvt",
"ndt", "hvt", "hvt", "hvt", "hvt", "hvt", "hvt", "ndt", "ndt",
"ndt", "ndt", "ndt", "mst", "ndt", "mst", "ndt", "ndt", "ndt",
"mst", "mst", "ndt", "hvt", "mst", "mst", "hvt", "mst", "hvt",
"ndt", "ndt", "hvt", "hvt", "mst", "mst", "hvt", "hvt", "mst",
"ndt", "mst", "hvt", "ndt", "mst", "mst", "mst", "mst", "ndt",
"mst", "hvt", "ndt", "hvt", "mst", "ndt", "hvt", "ndt", "ndt",
"ndt", "hvt", "ndt", "ndt", "mst", "mst", "mst", "hvt", "mst",
"mst", "ndt", "ndt", "hvt", "ndt", "ndt", "ndt", "hvt"), ASSOC = c("must",
"need_to", "need_to", "need_to", "need_to", "have_to", "need_to",
"need_to", "have_to", "must", "must", "need_to", "have_got_to",
"have_to", "must", "must", "have_to", "must", "have_to", "need_to",
"need_to", "need_to", "must", "have_to", "need_to", "have_to",
"need_to", "must", "need_to", "need_to", "have_got_to", "must",
"have_to", "need_to", "have_to", "must", "must", "must", "have_to",
"have_to", "have_to", "need_to", "must", "must", "need_to", "must",
"need_to", "have_to", "must", "have_to", "have_to", "must", "must",
"have_to", "need_to", "must", "need_to", "need_to", "need_to",
"have_to", "must", "need_to", "have_to", "have_to", "have_to",
"need_to", "need_to", "must", "need_to", "need_to", "have_to",
"have_to", "have_to", "have_to", "need_to", "must", "have_to",
"have_to", "need_to", "must", "have_to", "have_to", "must", "need_to",
"have_to", "have_to", "need_to", "must", "have_to", "must", "must",
"must", "must", "have_got_to", "need_to", "have_to", "need_to",
"need_to", "must", "have_to")), class = "data.frame", row.names = c(NA,
-100L))
df2[sapply(df2, is.character)] <- lapply(df2[sapply(df2, is.character)],
as.factor)
I have fitted the model like this:
multi_mo <- nnet::multinom(VERB ~ SOURC + SITN + PRSN + ANIM + ASSOC +
AGENT + CTGRY + CNTRY + CTGRY:CNTRY + SOURC:CTGRY + PRSN:CTGRY, data = df2, model=TRUE, maxit=1000)
Now I wanna see the marginal effects of variables SITN
and CNTRY
(following Sonderegger (2022): 176). Sonderegger uses ggeffects::ggeffect
but I also found marginaleffects::avg_prediction
while searching on the internet. However, I get different results from these two functions. For example see the predicted probabilities of dynamic
vs stative
for hvt
in panel 1 of each graph.
ggeffects::ggeffect
output:
library(ggeffects)
pred = ggeffect(multi_mo, terms = c("SITN", "CNTRY"))
p = ggplot(data=pred, aes(x=x, y=predicted)) +
facet_wrap(~ response.level) +
geom_point(aes(shape=group)) +
geom_line(aes(group = group, linetype = group))
print(p)
marginaleffects::avg_prediction
output
library(marginaleffects)
pred = avg_predictions (multi_mo, by = c("SITN", "CNTRY"), type = "probs")
p = ggplot(data=pred, aes(x=SITN, y=estimate)) +
facet_wrap(~ group) +
geom_point(aes(shape=CNTRY)) +
geom_line(aes(group = CNTRY, linetype = CNTRY))
print(p)
Am I comparing two different things here? If that is the case, what is the relevant function for ggeffects::ggeffect
in marginaleffects
package?