Showing posts with label R. Show all posts
Showing posts with label R. Show all posts

Wednesday, June 30, 2010

Multinomial logit in Stata and R III

Another set of translations between Stata and R - calculation of the most important kind of margins (see previous post), i.e., the marginal effects. This requires us to use yet another R package for out-of-the-box calculations of marginal effects, although I do also show how to do this by hand. Note that in both the R constructions, the standard errors are missing, and in a future post I will show how to calculate the standard errors analytically using the delta method and a nonparametric bootstrap.


Here I construct the marginal effects for the covariate "age" on the predicted probabilities of each of the three categories of the outcome. For the automatic construction of marginal effects in R, I make use of another package, VGAM and its vglm() and margeff() functions.
R

########################
# calculate the marginal effects
########################
library(VGAM)
library(gregmisc)
library(foreign)
sysdn <- read.dta("D:/programming/r/sandbox/sysdn1.dta", convert.factors=TRUE)
sysdn.multinomial <- vglm(insure~age+male+nonwhite+as.factor(site), data=sysdn, multinomial)
sysdn.multinomial
margeff.sysdn.multinomial <- margeff(sysdn.multinomial)
rowMeans(margeff.sysdn.multinomial["age",,])  # marginal effects
# manual construction of ME
p <- fitted(sysdn.multinomial)  # p_ij (row `i' = individual, column `j' = choice) here, a 644x3 matrix
avgmargeff.age1 <- mean(p[,1]*(coef(sysdn.multinomial)["age:1"] - p[,2]*coef(sysdn.multinomial)["age:2"]- p[,1]*coef(sysdn.multinomial)["age:1"]))
avgmargeff.age2 <- mean(p[,2]*(coef(sysdn.multinomial)["age:2"] - p[,2]*coef(sysdn.multinomial)["age:2"]- p[,1]*coef(sysdn.multinomial)["age:1"]))
# note that insure="Uninsure" is the base category and the coefficient vector is zero
avgmargeff.age3 <- mean(p[,3]*(0 - p[,2]*coef(sysdn.multinomial)["age:2"]- p[,1]*coef(sysdn.multinomial)["age:1"]))
avgmargeff.age1 
avgmargeff.age2
avgmargeff.age3
Stata
webuse sysdsn1, clear
mlogit insure age male nonwhite i.site, base(3)
margins, dydx(age) predict(outcome(1)) // average marginal effects
margins, dydx(age) predict(outcome(2))
margins, dydx(age) predict(outcome(3))
I think this back-and-forth between languages and packages makes a good case for learning more than one programming languages. Somethings are more easily done in one or the other - a valuable bit of flexibility when working on a large project. And to round it off, here is a very good explanation of what margins are, quoted from [R], Stata's base reference manual:
What we call margins of responses are also known as predictive margins, adjusted predictions, and recycled predictions. When applied to balanced data, margins of responses are also called estimated marginal means and least-squares means. A margin is a statistic based on a fitted model calculated over a dataset in which some of or all the covariates are fixed at values different from what they really are. For instance, after a linear regression fit on males and females, the marginal mean (margin of mean) for males is the predicted mean of the dependent variable, where every observation is treated as if it represents a male; thus those observations that in fact do represent males are included, as well as those observations that represent females. The marginal mean for female would be similarly obtained by treating all observations as if they represented females.

EDIT: Thank you to Gautam C. for pointing out what happens to HTML style comments if they are in the CSS of an HTML file.

Tuesday, June 29, 2010

Multinomial logit in Stata and R II

A couple more interesting translations:


  1. If you wanted to change the base outcome in the estimation, the multinom() function in nnet is no longer sufficient and we use the mlogit package. Note that this package is more in line with the motivation of multinomial logit as a choice model and thus expects the data in a form where for each choice occasion, all choices are enumerated as separate rows in the data.
    Stata
    webuse sysdsn1, clear
    mlogit insure age male nonwhite i.site, base(2)
    

    R
    library(effects)
    library(foreign)
    library(mlogit)
    # using the nnet package
    sysdn <- read.dta("D:/programming/r/sandbox/sysdn1.dta", convert.factors=TRUE)
    sysdn.multinom <- multinom(insure~age+male+nonwhite+as.factor(site), data=sysdn)
    summary(sysdn.multinom)
    # change the base category of the response
    sysdn$site <- as.factor(sysdn$site)
    sysdn$insure <- as.factor(sysdn$insure)
    sysdn.mldata <- mlogit.data(sysdn, varying=NULL, choice="insure", shape="wide")
    sysdn.mlogit <- mlogit(insure ~ 1|age+male+nonwhite+site, data=sysdn.mldata, reflevel="Indemnity") 
    sysdn.mlogit2 <- mlogit(insure ~ 1|age+male+nonwhite+site, data=sysdn.mldata, reflevel="Prepaid") 
    summary(sysdn.mlogit)
    summary(sysdn.mlogit2)
    

  2. You might want to predict (identified statistics) using the estimated model at various values of the covariates. In Stata 11, the new -margins- command offers a very general framework for doing this. I would strongly recommend reading [M] margins (the pdf documentation) for a very good description of what can be achieved using -margins-. Let us see is we can encode losslessly to R. I use the example I found here.
    R
    # calculate predicted probabilities
    mydata <- read.csv(url("http://www.ats.ucla.edu/stat/r/dae/mlogit.csv"))
    attach(mydata)
    names(mydata)
    mlogit.model<- mlogit(brand~1|female+age, data = mldata, reflevel="1")
    summary(mlogit.model)  
    newdata <- data.frame(cbind(age = rep(24:38, 2), female = c(rep(0, 15), rep(1, 15))))
    logit1 <- rep(0, 30)
    logit2 <- -11.774655 + 0.523814*newdata$female + 0.368206*newdata$age
    logit3 <- -22.721396 + 0.465941*newdata$female + 0.685908*newdata$age
    logits <- cbind(logit1, logit2, logit3)
    p.unscaled <- exp(logits)
    p <- cbind(newdata, (p.unscaled / rowSums(p.unscaled)))
    colnames(p) <- c("age", "female", "pred.1", "pred.2", "pred.3")
    p
    
    Stata
    insheet using http://www.ats.ucla.edu/stat/r/dae/mlogit.csv, comma clear
    mlogit brand female age, base(1)
    margins, predict(outcome(1)) at(age=(24(1)38) female=(0 1))
    margins, predict(outcome(2)) at(age=(24(1)38) female=(0 1))
    margins, predict(outcome(3)) at(age=(24(1)38) female=(0 1))
    
    The R code is more cumbersome, but it has the advantage of laying bare the structure of the construction, which is very valuable while learning about these models.

Multinomial logit in Stata and R

Here is how you can reproduce Stata's -mlogit- in R:
Stata


websuse sysdn1, clear
mlogit insure age male nowhite i.site

R

install.packages("foreign")
library(foreign)
sysdn <- read.dta("D:/programming/r/sandbox/sysdn1.dta", convert.factors=TRUE)
sysdn.multinom <- multinom(insure~age+male+nonwhite+ as.factor(site), data=sysdn)
summary(sysdn.multinom)

And as of today, I am using Revolution R's optimized 64-bit version of R which they have made free for academic use. The function multinom() in R comes from the pre-loaded nnet package. Syntax highlighting for R code here . Still looking for a brush/syntax highlighter for Stata. If anyone knows a good one, please do let me know.