R help - Apr 2016 - Predicting probabilities in ordinal probit analysis in R

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]]