########### # Functions to use the goodness-of-fit test of Bondell, # "Testing goodness-of-fit in logistic case-control studies", # Biometrika, 2007. ############ predict.prob<-function(x,alpha,beta) { exp(alpha+x%*%beta)/(1+exp(alpha+x%*%beta)) } dist2full <- function(dis) { n <- attr(dis, "Size") full <- matrix(0, n, n) full[lower.tri(full)] <- dis full + t(full) } C.matrix<-function(x,sigma) { e<-eigen(sigma) V<-e$vectors sigma.root<-V%*%diag(1/sqrt(e$values))%*%t(V) stand.x<-x%*%sigma.root temp.dist<-dist(stand.x) dnorm(dist2full(temp.dist), mean=0, sd=sqrt(2)) } ((sqrt(2)*sqrt(2*pi))^(1-p))* test.stat<-function(theta,x,y,C.quad.matrix) { n.0<-length(y[y==0]) n.1<-length(y[y==1]) rho.const<-n.1/n.0 alpha<-theta[1] beta<-theta[-1] resid.x<-y-predict.prob(x,alpha,beta) ((1+rho.const)/n.0)*t(resid.x)%*%C.quad.matrix%*%resid.x } dist.fit<-function(x,y) { fit.first<-glm(y~x,family=binomial,x=TRUE) theta<-as.vector(fit.first$coef) x0<-as.matrix(fit.first$x[,-1]) alpha<-theta[1] beta<-theta[-1] sigma.est<-cov(x0) C.quad.matrix<-C.matrix(x0,sigma.est) test.statistic<-test.stat(theta,x0,y,C.quad.matrix) return( list( test.statistic=test.statistic, theta=theta) ) } boot.dist<-function(x,y,alpha,beta) { quad.test.stat.star<-NULL B<-2000 n.0<-length(y[y==0]) n.1<-length(y[y==1]) for (j in 1:B) { num.vec<-1:length(x[,1]) x.0.star.temp<-sample(num.vec, n.0, replace = TRUE, prob = (n.0*(1+exp(alpha+x%*%beta)))^(-1)) x.1.star.temp<-sample(num.vec, n.1, replace = TRUE, prob = exp(alpha+x%*%beta)*(n.0*(1+exp(alpha+x%*%beta)))^(-1)) x.0.star<-as.matrix(x[x.0.star.temp,]) x.1.star<-as.matrix(x[x.1.star.temp,]) x.star<-as.matrix(rbind(x.0.star,x.1.star)) y.star<-c(rep(0,n.0),rep(1,n.1)) quad.test.stat.star[j]<-dist.fit(x.star,y.star)$test.statistic } sort(quad.test.stat.star) } min.dist.test<-function(x,y) { fitted.stat<-dist.fit(x,y) test.statistic<-as.numeric(fitted.stat$test.statistic) theta<-fitted.stat$theta alpha<-theta[1] beta<-theta[-1] null.dist<-boot.dist(x,y,alpha,beta) p.value.test<-sum(null.dist>test.statistic)/length(null.dist) return( list( test.statistic=test.statistic, null.dist=null.dist, theta=theta, p.value.test=p.value.test) ) }