[R] Standard Error for difference in predicted probabilities

[Darin A. England]

I think you can use the bootstrap to obtain the std error. My
attempt for your problem and data is below. I would be interested if
anyone can point out a problem with this approach.
Darin

y=rbinom(100,1,.4)
x1=rnorm(100, 3, 2)
x2=rbinom(100, 1, .7) 

diff <- vector(mode="numeric", length=200)
for (i in 1:200) {
    idx <- sample(1:length(y), length(y), replace=T)
    mod=glm(y[idx] ~ x1[idx] + x2[idx], family=binomial)
    pred=predict(mod, newdata=data.frame(cbind(x1=rep(mean(x1[idx]),
        100), x2=x2[idx])), type="response")
    diff[i]=unique(pred)[1]-unique(pred)[2]
}
sd(diff)


On Fri, Sep 24, 2010 at 09:17:29AM -0400, Andrew Miles wrote:
> Is there a way to estimate the standard error for the difference in  
> predicted probabilities obtained from a logistic regression model?
> 
> For example, this code gives the difference for the predicted  
> probability of when x2==1 vs. when x2==0, holding x1 constant at its  
> mean:
> 
> y=rbinom(100,1,.4)
> x1=rnorm(100, 3, 2)
> x2=rbinom(100, 1, .7)
> mod=glm(y ~ x1 + x2, family=binomial)
> pred=predict(mod, newdata=data.frame(cbind(x1=rep(mean(x1), 100),  
> x2)), type="response")
> diff=unique(pred)[1]-unique(pred)[2]
> diff
> 
> I know that predict() will output SE's for each predicted value, but  
> how do I generate a SE for the difference in those predicted values?
> 
> Thanks in advance!
> 
> Andrew Miles
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.