[R] Predicting probabilities in ordinal probit analysis in R

[Faradj Koliev]

Dear all, 

I have two questions that are almost completely related to how to do things in R.

I am running an ordinal probit regression analysis in R. The dependent variable has three levels (0=no action; 1=warning; 2=sanction).

I use the lrm command in the rms package:

print( res1<- lrm(Y ~ x1+x2+x3+x4+x5+x6, y=TRUE, x=TRUE, data=mydata))
I simply couldn't make any sense of the information generated my ?predict.lrm. What I want to do is to calculate the marginal effects of all explanatory variables for each level of the dependent variable. In Stata, this is very simple: mfx compute, predict (outcome(#0)); mfx compute, predict (outcome(#2)) and mfx compute, predict (outcome(#3)).

So my first question is: how do I generate marginal effects for each outcome in R? 

The second question is related to interaction effects, which I need to include in the same model:

print( res1<- lrm(Y ~ x1+x2+x3+x4+x5+x6+x5*x6, y=TRUE, x=TRUE, data=mydata))
If I knew the answer to the first question, I would have ran marginal effects with the interaction term included. Then, I would have plotted the predicted values of the interaction term.

So the second question is: how do I plot the effects (predicted values) of variables in the interaction term?

Many thanks!

Small sample from my dataset (only one country)

dput(mydatasample)

structure(list(year = 1989:2014, country = structure(c(1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "Canada", class = "factor"), 
    id = structure(1:26, .Label = c("CAN 1989", "CAN 1990", "CAN 1991", 
    "CAN 1992", "CAN 1993", "CAN 1994", "CAN 1995", "CAN 1996", 
    "CAN 1997", "CAN 1998", "CAN 1999", "CAN 2000", "CAN 2001", 
    "CAN 2002", "CAN 2003", "CAN 2004", "CAN 2005", "CAN 2006", 
    "CAN 2007", "CAN 2008", "CAN 2009", "CAN 2010", "CAN 2011", 
    "CAN 2012", "CAN 2013", "CAN 2014"), class = "factor"), stage1 = c(1L, 
    1L, 0L, 0L, 0L, 0L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 0L, 0L, 0L, 
    0L, 2L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L), x1 = c(1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L), x2 = c(1L, 2L, 1L, 2L, 2L, 
    2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 
    2L, 1L, 2L, 2L, 2L, 2L), x3 = c(9L, 9L, 9L, 9L, 9L, 9L, 9L, 
    9L, 9L, 9L, 9L, 8L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
    10L, 10L, 10L, 10L, 10L, 10L), x4 = c(31L, 31L, 31L, 31L, 
    31L, 30L, 30L, 30L, 31L, 30L, 29L, 30L, 28L, 28L, 28L, 27L, 
    29L, 29L, 29L, 28L, 25L, 24L, 23L, NA, NA, NA), x5 = structure(1:26, .Label = c("17,12528685", 
    "17,14022279", "17,15382785", "17,16610202", "17,17704534", 
    "17,18665779", "17,19493938", "17,20571103", "17,21628118", 
    "17,22493732", "17,23321101", "17,242041", "17,25213621", 
    "17,26110753", "17,27106985", "17,2810902", "17,29094924", 
    "17,29891768", "17,30861622", "17,31943819", "17,33088659", 
    "17,34202619", "17,35190237", "17,36381421", "17,37537139", 
    "17,38618117"), class = "factor"), x5.1 = c(0L, 0L, 0L, 0L, 
    1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 
    0L, 0L, 0L, 0L, 1L, 1L, 0L)), .Names = c("year", "country", 
"id", "stage1", "x1", "x2", "x3", "x4", "x5", "x5.1"), class = "data.frame", row.names = c(NA, 
-26L))




	[[alternative HTML version deleted]]