Edit 2: I changed the link (earlier was the wrong link)
Edit 1: My dataframe for this question is like:
| 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 |
continued:
| 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 |
I want to make effects plot of my polr model as they do in this example: https://data.library.virginia.edu/visualizing-the-effects-of-proportional-odds-logistic-regression/
I tried and got error which I don't understand why and what I do wrong. My codes:
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)
The above codes works and gives the outputs.
The below code does not work and give error:
plot(Effect(focal.predictors = c("country","work less"), mod = for.plot, xlevels = list(age = 15:65)), rug = FALSE)
this gives the following error:
Error in \contrasts<-(\tmp, value = contr.funs[1 + \isOF[nn]]):
contrasts can be applied only to factors with 2 or more levels
Does anyone know how I can solve this error?
structure of my data:
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 ...
New try - with as.factor:
running the regression works fine, and it gives the summary and anova output:
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)
running this give a new error:
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
New try again - making the variables to factor. The structur of the data:
$ 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 ...
The structure of the data in the example is: 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 ...
But it still doesn't work :/
Okay, now I figured out what was wrong. I need to recode the variables using as.factor to change them in the dataframe instead of using the as.factor inside the polr regression.
So, I code it as:
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)
and then running the regression with:
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)
Gives me a plot very much similar to this:
country * age * gender effect plot
