#The function R2.adj computes the two adjusted coefficients of determination 
#for logistic regression described in Liao and McGee (2003) Adjusted Coefficients 
#of determination for logistic regression. The American Statistician, August issue

#########################################################
#Example:
#  y = rbinom(100, 1, .5)     #Note that one subject must be on a separate row
#  x = rnorm(500)
#  dim(x) = c(100, 5)
#  R2.adj(y~x) # computes the adjusted coefficients of determination,
               # the first is for the entropy loss and the second for the square error loss
#
#  R2(y~x)     #computes the unadjusted coefficients of determination
###########################################################


logistic.fitted = function(y,x) 
{
   obj = glm(y~x, family=binomial);
   pi = obj$fitted.values;

   lower.limit = .00001;
   upper.limit = .99999;
   pi[pi<lower.limit] = lower.limit;
   pi[pi>upper.limit] = upper.limit;

   pi;
}


log.likelihood.null = function(y)
{
   fitted = mean(y);

   value1 = -mean( y*log(fitted) + (1-y)*log(1-fitted) );
   value2 = mean( (y-fitted)^2 )

   c(value1, value2);
}

log.likelihood.null.true = function(y, pi)
{
   value1 = -mean( y*log(pi) + (1-y)*log(1-pi) );
   value2 = mean( (y-pi)^2 )

   c(value1, value2);
}

log.likelihood = function(y, x)
{
   fitted = logistic.fitted(y, x);

   value1 = -mean( y*log(fitted) + (1-y)*log(1-fitted) );
   value2 = mean( (y-fitted)^2 )

   c(value1, value2);
}

log.likelihood.true = function(y, pi)
{
   value1 = -mean( y*log(pi) + (1-y)*log(1-pi) );
   value2 = mean( (y-pi)^2 )

   c(value1, value2);
}

#####################################################################
bias = function(pi, x)
{
   ave = numeric(2);
   n.simu = 500;
   for(i in 1:n.simu)
   {
       y = rbinom(length(pi), 1, pi);
       ave = ave + log.likelihood(y, x) - log.likelihood.true(y, pi);
   }

   ave/n.simu;
}

bias.null = function(pi, n)   #pi is a scalar
{
   ave = numeric(2);
   n.simu = 500;
   for(i in 1:n.simu)
   {
       y = rbinom(n, 1, pi);
       ave = ave + log.likelihood.null(y) - log.likelihood.null.true(y, pi);
   }

   ave/n.simu;
}

###################################################################
R2.adj = function(formula)
{
   obj = glm(formula, family=binomial, x=TRUE);
   y = obj$y;
   x = obj$x[,-1]; #without the column of 1

   fitted = logistic.fitted(y, x);
   ll.model = log.likelihood(y, x) - bias(fitted, x);

   ll.null = log.likelihood.null(y) - bias.null(mean(y), length(y));

   value = 1 - ll.model/ll.null;
   value[value<0] = 0;

   value;
 }
###################################################################
R2 = function(formula)
{
   obj = glm(formula, family=binomial, x=TRUE);
   y = obj$y;
   x = obj$x[,-1]; #without the column of 1

   1- log.likelihood(y, x)/log.likelihood.null(y);
}