wlogreg <-
function(x0,y,initwml=FALSE,const=0.5,kmax=1e3,maxhalf=10)
{
# Computation of the estimator of Bianco and Yohai (1996) in logistic regression
# -------------
# This is a slightly modified version of code due to
# Christophe Croux, Gentiane Haesbroeck, and Kristel Joossens
# (Here initwml defaults to F
#
# This program computes the estimator of Bianco and Yohai in
# logistic regression. By default, an intercept term is included
# and p parameters are estimated.
#
# For more details we refer to
# Croux, C., and Haesbroeck, G. (2003), ``Implementing the Bianco and Yohai
# estimator for Logistic Regression'',
# Computational Statistics and Data Analysis, 44, 273-295
#
#Input:
#-------
# x0= n x (p-1) matrix containing the explanatory variables;
# y= n-vector containing binomial response (0 or 1);
#
# initwml= logical value for selecting one of the two possible methods for computing
# the initial value of the optimization process. If initwml=T (default), a
# weighted ML estimator is computed with weights derived from the MCD estimator
# computed on the explanatory variables. If initwml=F, a classical ML fit is perfomed.
# When the explanatory variables contain binary observations, it is recommended
# to set initwml to F or to modify the code of the algorithm to compute the weights
# only on the continuous variables.
# const= tuning constant used in the computation of the estimator (default=0.5);
# kmax= maximum number of iterations before convergence (default=1000);
# maxhalf= max number of step-halving (default=10).
#
# Example:
# x0=matrix(rnorm(100,1))
# y0=numeric(runif(100)>0.5)
# BYlogreg(x0,y)
#
#Output:
#--------
# list with
# 1st component: T or F if convergence achieved or not
# 2nd component: value of the objective function at the minimum
# p next components: estimates for the parameters.
# p last components: standard errors of the parameters (if first component is T)
library(MASS)
x0=as.matrix(x0)
# n=nrow(x0)
p=ncol(x0)+1
p0=p-1
#Smallest value of the scale parameter before implosion
sigmamin=1e-4
# eliminate any rows with missing values
zz=elimna(cbind(x,y))
x=as.matrix(zz[,1:p0])
y=zz[,p]
n=nrow(x)
print(n)
# x=as.matrix(cbind(rep(1,n),x0))
x=as.matrix(cbind(rep(1,n),x))
print(rep(1,n))
y=as.numeric(y)
# Computation of the initial value of the optimization process
if (initwml==T)
{
hp=floor(n*(1-0.25))+1
mcdx=cov.mcd(x0, quantile.used =hp,method="mcd")
rdx=sqrt(mahalanobis(x0,center=mcdx$center,cov=mcdx$cov))
vc=sqrt(qchisq(0.975,p-1))
wrd=(rdx<=vc)
gstart=glm(y~x0,family=binomial,subset=wrd)$coef
}
else {gstart=glm(y~x0,family=binomial)$coef}
sigmastart=1/sqrt(sum(gstart^2))
xistart=gstart*sigmastart
stscores=x %*% xistart
sigma1=sigmastart
#Initial value for the objective function
oldobj=mean(phiBY3(stscores/sigmastart,y,const))
kstep=jhalf=1
while ((kstep < kmax) & (jhalf<maxhalf)){
unisig <- function(sigma)
{ mean(phiBY3(stscores/sigma,y,const))}
optimsig=nlminb(sigma1,unisig,lower=0)
sigma1=optimsig$par
if (sigma1<sigmamin) {print("Explosion");kstep=kmax
} else {
gamma1=xistart/sigma1
scores=stscores/sigma1
newobj=mean(phiBY3(scores,y,const))
oldobj=newobj
gradBY3=colMeans((derphiBY3(scores,y,const)%*%matrix(1,ncol=p))*x)
h=-gradBY3+((gradBY3 %*% xistart) *xistart)
finalstep=h/sqrt(sum(h^2))
xi1=xistart+finalstep
xi1=xi1/(sum(xi1^2))
scores1=(x%*%xi1)/sigma1
newobj=mean(phiBY3(scores1,y,const))
####stephalving
hstep=jhalf=1
while ((jhalf <=maxhalf) & (newobj>oldobj)){
hstep=hstep/2
xi1=xistart+finalstep*hstep
xi1=xi1/sqrt(sum(xi1^2))
scores1=x%*%xi1/sigma1
newobj=mean(phiBY3(scores1,y,const))
jhalf=jhalf+1
}
CONV=F
if ((jhalf==maxhalf+1) & (newobj>oldobj)) {CONV=T
} else {
jhalf=1
xistart=xi1
oldobj=newobj
stscores=x%*% xi1
kstep=kstep+1
}
}
}
if (kstep == kmax) {
CONV=F # print("No convergence")
result=list(convergence=FALSE,objective=0,coef=t(rep(NA,p)))
} else {
gammaest=xistart/sigma1
stander=sterby3(x0,y,const,gammaest)
result=list(convergence=CONV,coef=t(gammaest),sterror=stander)
}
return(result)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.