어떻게 LM 후 R() STATA에서

Question

에서 "에서 여백"STATA의 복제합니다. 여백 여백은 전체 여백 - 일정한 것이없는 여백 -을보고합니다. 우리 모델 이 물류 적이기 때문에 예측 된 확률의 평균값이보고됩니다. at() 옵션은 지정된 값으로 하나 이상의 공변량을 으로 수정하고 변수 및 연속 변수와 함께 사용할 수 있습니다. 마진이 나이 = 40을 설정, 데이터 모두에 대한 응답을 평균 이상으로 할 다음 따라서, 경우에 당신은

margins, at(age=40) 

를 입력했습니다.

어떤 패키지가 유용 할 수 있습니까? 하위 데이터에 대한 예측 값의 평균을 찾기 위해 이미 시도했지만 시퀀스 (예 : 여백 = 40 (1) 50)에서는 작동하지 않습니다.

Answer 1

당신은 STATA의 margins, at수단 또는 대표자 점 (this 및 설명서 참조)에서 평가 한계 효과는 단순히 이해한다 R.

에 한계 효과를 얻을 수있는 여러 가지 방법이 있습니다.

나는 그것이 가장 유사대로 사용하는지 가장이 솔루션을 좋아하고 있다고 생각 :

library(devtools) 
install_github("leeper/margins") 

출처 : https://github.com/leeper/margins

마진 포트 STATA의에 대한 노력 (이다 폐색 소스) 마진 명령 을 모델 객체에 포함 된 공변량의 경계 효과 ( "부분 효과") (예 : 클래스 "lm"및 "glm"의 )를 계산하기위한 S3 일반 방법으로 사용합니다. 새 "여백"클래스의 플롯 방법 은 marginsplot 명령을 추가로 포트합니다.

 cyl  hp  wt 
0.03814 -0.04632 -3.11981 

library(margins) 
x <- lm(mpg ~ cyl * hp + wt, data = mtcars) 
(m <- margins(x)) 

이 패키지 또한 prediction 명령 ( ?prediction)를 참조하십시오.

I. erer (패키지)

maBina() command 

http://cran.r-project.org/web/packages/erer/erer.pdf

II : 그에서

잡담은 여기에 몇 가지 다른 내가 컴파일 한 솔루션입니다. http://www.r-bloggers.com/probitlogit-marginal-effects-in-r/

III :

mfxboot <- function(modform,dist,data,boot=1000,digits=3){ 
    x <- glm(modform, family=binomial(link=dist),data) 
    # get marginal effects 
    pdf <- ifelse(dist=="probit", 
       mean(dnorm(predict(x, type = "link"))), 
       mean(dlogis(predict(x, type = "link")))) 
    marginal.effects <- pdf*coef(x) 
    # start bootstrap 
    bootvals <- matrix(rep(NA,boot*length(coef(x))), nrow=boot) 
    set.seed(1111) 
    for(i in 1:boot){ 
    samp1 <- data[sample(1:dim(data)[1],replace=T,dim(data)[1]),] 
    x1 <- glm(modform, family=binomial(link=dist),samp1) 
    pdf1 <- ifelse(dist=="probit", 
        mean(dnorm(predict(x, type = "link"))), 
        mean(dlogis(predict(x, type = "link")))) 
    bootvals[i,] <- pdf1*coef(x1) 
    } 
    res <- cbind(marginal.effects,apply(bootvals,2,sd),marginal.effects/apply(bootvals,2,sd)) 
    if(names(x$coefficients[1])=="(Intercept)"){ 
    res1 <- res[2:nrow(res),] 
    res2 <- matrix(as.numeric(sprintf(paste("%.",paste(digits,"f",sep=""),sep=""),res1)),nrow=dim(res1)[1]) 
    rownames(res2) <- rownames(res1) 
    } else { 
    res2 <- matrix(as.numeric(sprintf(paste("%.",paste(digits,"f",sep=""),sep="")),nrow=dim(res)[1])) 
    rownames(res2) <- rownames(res) 
    } 
    colnames(res2) <- c("marginal.effect","standard.error","z.ratio") 
    return(res2)} 

소스를 mfxboot.출처 : 지금 R probit regression marginal effects

x1 = rbinom(100,1,.5) 
x2 = rbinom(100,1,.3) 
x3 = rbinom(100,1,.9) 
ystar = -.5 + x1 + x2 - x3 + rnorm(100) 
y = ifelse(ystar>0,1,0) 
probit = glm(y~x1 + x2 + x3, family=binomial(link='probit')) 
xbar <- as.matrix(mean(cbind(1,ttt[1:3]))) 

그래픽, 즉, X1, X2 및 X3의 한계 효과

library(arm) 
curve(invlogit(1.6*(probit$coef[1] + probit$coef[2]*x + probit$coef[3]*xbar[3] + probit$coef[4]*xbar[4]))) #x1 
curve(invlogit(1.6*(probit$coef[1] + probit$coef[2]*xbar[2] + probit$coef[3]*x + probit$coef[4]*xbar[4]))) #x2 
curve(invlogit(1.6*(probit$coef[1] + probit$coef[2]*xbar[2] + probit$coef[3]*xbar[3] + probit$coef[4]*x))) #x3 

library(AER) 
data(SwissLabor) 
mfx1 <- mfxboot(participation ~ . + I(age^2),"probit",SwissLabor) 
mfx2 <- mfxboot(participation ~ . + I(age^2),"logit",SwissLabor) 
mfx3 <- mfxboot(participation ~ . + I(age^2),"probit",SwissLabor,boot=100,digits=4) 

mfxdat <- data.frame(cbind(rownames(mfx1),mfx1)) 
mfxdat$me <- as.numeric(as.character(mfxdat$marginal.effect)) 
mfxdat$se <- as.numeric(as.character(mfxdat$standard.error)) 


# coefplot 
library(ggplot2) 
ggplot(mfxdat, aes(V1, marginal.effect,ymin = me - 2*se,ymax= me + 2*se)) + 
    scale_x_discrete('Variable') + 
    scale_y_continuous('Marginal Effect',limits=c(-0.5,1)) + 
    theme_bw() + 
    geom_errorbar(aes(x = V1, y = me),size=.3,width=.2) + 
    geom_point(aes(x = V1, y = me)) + 
    geom_hline(yintercept=0) + 
    coord_flip() + 
    opts(title="Marginal Effects with 95% Confidence Intervals")