绘制 polr 模型的效果时如何修复错误？

Question

编辑 2：我更改了链接（之前是错误的链接）
编辑 1：我对这个问题的数据框是这样的：

gender	Work less	happy	lifestatisfied	country	Work much
2	0	7	8	GB	1
1	1	8	8	SE	0
1	0	7	9	DK	1
1	0	6	9	DE	1
1	NA	7	5	NO	NA

继续：

health	education	income	age	marital status
3	3	Na	61	NA
4	2	2	30	NA
1	3	4	39	6
5	7	5	52	4
4	1	5	17	5

• country 是字符（即国家名称）我有 5 个国家
• 性别是虚拟的 1 或 2
• 年龄是受访者的年龄，例如 35、47 等。
• 收入按比例缩放，为 1 到 10
• educ（教育）是1到7
• 健康等级从 1 到 5
• 工作少是虚拟的，即 1 或 0
• work much 是虚拟的，即 1 或 0
• 婚姻状况从 1 到 6

我想像他们在这个例子中那样制作我的 polr 模型的效果图：https://data.library.virginia.edu/visualizing-the-effects-of-proportional-odds-logistic-regression/

我尝试并得到了错误，我不明白为什么以及我做错了什么。 我的代码：

for.plot <- polr(factor(as.ordered(lifesatisfaction)) ~ country*(gender + age + income + educ + health + work less + work much), data = surveywave5, method = "logistic", Hess = TRUE)

summary(for.plot)

Anova(for.plot)

上面的代码有效并给出了输出。

下面的代码不起作用并给出错误：

plot(Effect(focal.predictors = c("country","work less"), mod = for.plot, xlevels = list(age = 15:65)), rug = FALSE)

这给出了以下错误：

Error in \contrasts<-(\tmp, value = contr.funs[1 + \isOF[nn]]):

contrasts can be applied only to factors with 2 or more levels

有谁知道我该如何解决这个错误？

我的数据结构：

str(surveywave5)
tibble [7,936 x 25] (S3: tbl_df/tbl/data.frame)
 $ cntry            : chr+lbl [1:7936] DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE, DE...
   ..@ label       : chr "Country"
   ..@ format.stata: chr "%7s"
   ..@ labels      : Named num [1:27] 0 0 0 0 0 0 0 0 0 0 ...
   .. ..- attr(*, "names")= chr [1:27] "Denmark" "United Kingdom" "Norway" "Germany" ...
 $ gender           : dbl+lbl [1:7936] 2, 1, 1, 1, 1, 1, 2, 1, 2, 1, 2, 2, 2, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 2, 2, 1, 2, 2, 2, 2, 1, 2, 1, 2, 1, 1, 2, 1, 1, 2, 2, 1, 2, 2, ...
   ..@ label       : chr "Gender"
   ..@ format.stata: chr "%6.0g"
   ..@ labels      : Named num [1:3] 1 2 9
   .. ..- attr(*, "names")= chr [1:3] "Male" "Female" "No answer"
 $ marital status   : num [1:7936] NA NA NA 6 6 NA NA NA NA 6 ...
 $ age             : dbl+lbl [1:7936] 61, 30, 39, 52, 17, 50, 31, 49, 50, 39, 45, 42, 62, 27, 52, 19, 17, 20, 62, 34, 18, 59, 17, 46, 37, 59, 49, 21, 27, 40, 47, 52, 21, 46, 38...
   ..@ label       : chr "Age of respondent, calculated"
   ..@ format.stata: chr "%8.0g"
   ..@ labels      : Named num 999
   .. ..- attr(*, "names")= chr "Not available"
 $ income           : num [1:7936] NA 2 4 5 5 10 7 7 6 3 ...
 $ education        : num [1:7936] 3 2 3 7 1 7 3 3 3 5 ...
 $ lifesatisfaction : num [1:7936] 8 8 9 9 9 8 10 9 10 6 ...
 $ happy            : num [1:7936] 9 8 8 8 10 7 10 10 9 6 ...
 $ health           : num [1:7936] 3 4 1 3 4 5 5 4 4 3 ...
 $ work less        : num [1:7936] 0 1 0 0 NA 0 0 0 1 0 ...
 $ work much        : num [1:7936] 1 0 1 1 NA 0 1 1 0 1 ...
 $ marital status   : num [1:7936] 1 1 1 NA NA 1 1 1 1 NA ...

新尝试 - 使用 as.factor：

运行回归工作正常，它给出了摘要和方差分析输出：

for.plot <- polr(as.factor(as.ordered(lifesatisfaction)) ~ as.factor(country)*(gender + age + income + education + health + as.factor(work less) + work much), data = surveywave5, method = "logistic", Hess = TRUE)

运行它会产生一个新错误：

plot(Effect(focal.predictors = c("country","work less"), mod = for.plot, xlevels = list(age = 15:65)), rug = FALSE)

Error in Effect(focal.predictors = c("country", "work less"), mod = for.plot, : 
model formula should not contain calls to
factor(), as.factor(), ordered(), as.ordered(), as.numeric(), or as.integer();
see 'Warnings and Limitations' in ?Effect

再试一次 - 使变量成为因素。数据结构：

 $ country.new           : Factor w/ 5 levels "DE","DK","GB",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ lifesatisfaction.new  : Ord.factor w/ 11 levels "0"<"1"<"2"<"3"<..: 9 9 10 10 10 9 11 10 11 7 ...        
 $ gender.new          : Factor w/ 2 levels "1","2": 2 1 1 1 1 1 2 1 2 1 ...
 $ income.new      : Factor w/ 10 levels "1","2","3","4",..: NA 2 4 5 5 10 7 7 6 3 ...
 $ eisced          : Factor w/ 7 levels "1","2","3","4",..: 3 2 3 7 1 7 3 3 3 5 ...
 $ health.new        : Factor w/ 5 levels "1","2","3","4",..: 3 4 1 3 4 5 5 4 4 3 ...
 $ work.less.new  : Factor w/ 2 levels "0","1": 1 2 1 1 NA 1 1 1 2 1 ...
 $ work much.new  : Factor w/ 2 levels "0","1": 2 1 2 2 NA 1 2 2 1 2 ...  

示例中的数据结构为： Visualizing the Effects of Proportional-Odds Logistic Regression

str(WVS)
'data.frame':   5381 obs. of  6 variables:
 $ poverty : Ord.factor w/ 3 levels "Too Little"<"About Right"<..: 1 2 1 3 1 2 3 1 1 1 ...
 $ religion: Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
 $ degree  : Factor w/ 2 levels "no","yes": 1 1 1 2 2 1 1 1 1 1 ...
 $ country : Factor w/ 4 levels "Australia","Norway",..: 4 4 4 4 4 4 4 4 4 4 ...
 $ age     : int  44 40 36 25 39 80 48 32 74 30 ...
 $ gender  : Factor w/ 2 levels "female","male": 2 1 1 1 2 1 1 2 1 2 ...

但它仍然不起作用：/

Answer 1

好的，现在我知道出了什么问题。我需要使用 as.factor 重新编码变量以在数据框中更改它们，而不是在 polr 回归中使用 as.factor 。
所以，我把它编码为：

surveywave5$country = as.factor(surveywave5$country)
surveywave5$lifesatisfied = as.factor(as.ordered((surveywave5$lifesatisfied)))
surveywave5$gender = as.factor(surveywave5$gender)
surveywave5$age = as.integer(surveywave5$age)
surveywave5$income = as.factor(surveywave5$income)
surveywave5$work less = as.factor(surveywave5$work less)  

然后运行回归：

for.plot <- polr(lifesatisfaction ~ country*(gender + age + income + health + work less), data = surveywave5, method = "logistic", Hess = TRUE)

plot(Effect(focal.predictors = c("country","age","gender"), mod = for.plot, xlevels = list(age = 15:65), latent = TRUE), rug = FALSE)

给我一个与此非常相似的情节：
country * age * gender effect plot