Мне не понятно, что ты пытаешься сделать. glm(..., family = binomial(link = "logit"))
соответствует простой логистической модели, предполагающей бинарный результат (например, успех против неудачи, следовательно, базовое биномиальное распределение).
Если вы хотите установить модель с полиномом , вы можете использовать nnet::multinom
library(nnet)
fit <- multinom(Species ~ ., data = iris)
#Call:
#multinom(formula = Species ~ ., data = iris)
#
#Coefficients:
# (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
#versicolor 18.69037 -5.458424 -8.707401 14.24477 -3.097684
#virginica -23.83628 -7.923634 -15.370769 23.65978 15.135301
#
#Residual Deviance: 11.89973
#AIC: 31.89973