Biost 536 Oct 4 lecture in R

Oct4.R

##biost 536, University of Washington
##this code generally corresponds to the stata code in the October 4th lecture
d <- expand.grid(case=0:1, xray=0:1,yob=0:1)
d$freq <- c(2388,2220,241,409,2936,2774,361,523)


##model with no interactions
mylogit <- glm(case~xray + yob, data=d,
	family=binomial(link="logit"),
	weights=d$freq)

confint.default(mylogit) ##normal confidence intervals--identical to stata
confint(mylogit) ##will provide profile likelihood confidence intervals 

#convenience function to get exponentiated output similar to stata
#you'll still need summary(x) for some information (degrees of freedom, etc)
#I couldn't figure out how to reproduce the standard error of the 
#exponentiated coefficient that stata reports


##x is the glm model object
e <- function(x, decimals=2, p.decimals=3){
	cbind(
		OR = round(exp(coef(x)),decimals),
		z=round(coef(summary(x))[,'z value'],decimals),
		'Pr(>|z|)'=round(coef(summary(x))[,'Pr(>|z|)'],
			p.decimals),
		round(exp(confint.default(x)),decimals)
	)
}

e(mylogit)


#with interactions
mylogit2 <- glm(case~xray*yob, data=d,
	family=binomial(link="logit"),
	weights=d$freq)


e(mylogit2,5)

###OR for older group can be calculated by hand as
1.8255*.839


###or using lincom equivalents


require(gmodels)
###confidence intervals are not exactly the same as stata, but close 
gmodels.lincom <- estimable(mylogit2, c("xray"=1, "xray:yob"=1),conf=.95)
exp(gmodels.lincom)[,c(1,6,7)]

###identical to stata--find confidence intervals using the typical method of
## exp(estimate+1.96*SE)
exp(c(gmodels.lincom$Estimate, gmodels.lincom$`Std. Error`) %*% 
	c(1, -qnorm(.975)))
exp(c(gmodels.lincom$Estimate, gmodels.lincom$`Std. Error`) %*% 
	c(1, qnorm(.975)))


##an alternate method using a different package, identical to stata

require(survey)
temp <- svycontrast(mylogit2,quote(xray + `xray:yob`))
#note that default print methods of temp don't reflect str(temp)
survey.lincom <- c(temp, "se"=sqrt(attr(temp,"var")))

exp(survey.lincom %*% c(1, -qnorm(.975)))
exp(survey.lincom %*% c(1, qnorm(.975)))

##also see http://www.mail-archive.com/r-help@r-project.org/msg40301.html

###convenience wrapper function
##no guarantee this always works...

lincom <- function(x, expr,decimals=2, p.decimals=3){
	temp <- svycontrast(x,substitute(expr))
	survey.lincom <- as.vector(c(temp, "se"=sqrt(attr(temp,"var"))))
	
	c(
		Estimate=round(exp(survey.lincom[1]),decimals),
		z=round(survey.lincom[1]/survey.lincom[2],decimals),
		"Pr(>|z|)"=round(pnorm(
			abs(survey.lincom[1]/survey.lincom[2]),
			lower.tail=FALSE)*2,p.decimals),
		"2.5 %"=round(exp(survey.lincom %*% 
			c(1, -qnorm(.975))),decimals),
		"97.5 %"=round(exp(survey.lincom %*% 
			c(1, qnorm(.975))),decimals)
	)
}
lincom(x,xray+`xray:yob`,3)

##using reparameterization

d$yob2 <- as.integer(!d$yob)

mylogit3 <- glm(case~xray*yob2, data=d,
	family=binomial(link="logit"),
	weights=d$freq)

e(mylogit3,3)