Handcoding a Logit Model Below is an example of how to handcode a logit model. set.seed(123) n=100 N=1+rpois(n,5) X1=runif(n) X2=rexp(n) s=X2-X1-2 p=exp(s)/(1+exp(s)) vY=NULL for(i in 1:n){ Y=rbinom(1,prob=p[i],size=N[i]) vY=c(vY,Y) } db=data.frame(Y=vY,N=N,X1,X2) head(db,4) ## Y N X1 X2 ## 1 3 5 0.5999890 1.80312932 ## 2 1 8 0.3328235 0.03005945 ## 3 0 5 0.4886130 1.30344055 ## 4 0 9 0.9544738 0.19979224 vY=vX1=vX2=vN=NULL for(i in 1:n){ vY=c(vY,c(rep(0,db$N[i]-db$Y[i]),rep(1,db$Y[i]))) vX1=c(vX1,rep(db$X1[i],db$N[i])) vX2=c(vX2,rep(db$X2[i],db$N[i])) } largedb=data.frame(Z=vY,X1=vX1,X2=vX2) head(largedb,16) ## Z X1 X2 ## 1 0 0.5999890 1.80312932 ## 2 0 0.5999890 1.80312932 ## 3 1 0.5999890 1.80312932 ## 4 1 0.5999890 1.80312932 ## 5 1 0.5999890 1.80312932 ## 6 0 0.3328235 0.03005945 ## 7 0 0.3328235 0.03005945 ## 8 0 0.3328235 0.03005945 ## 9 0 0.3328235 0.03005945 ## 10 0 0.3328235 0.03005945 ## 11 0 0.3328235 0.03005945 ## 12 0 0.3328235 0.03005945 ## 13 1 0.3328235 0.03005945 ## 14 0 0.4886130 1.30344055 ## 15 0 0.4886130 1.30344055 ## 16 0 0.4886130 1.30344055 hours <- c(0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00, 3.25, 3.50, 4.00, 4.25, 4.50, 4.75, 5.00, 5.50) pass <- c( 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1) glm(pass ~ hours, family = binomial) ## ## Call: glm(formula = pass ~ hours, family = binomial) ## ## Coefficients: ## (Intercept) hours ## -4.078 1.505 ## ## Degrees of Freedom: 19 Total (i.e. Null); 18 Residual ## Null Deviance: 27.73 ## Residual Deviance: 16.06 AIC: 20.06 # for two hours 1 / ( 1 + exp(-(-4.0777+1.5046*2) ) ) ## [1] 0.2556884 download.file('http://www.ats.ucla.edu/stat/sas/dae/binary.sas7bdat', destfile = "binary.sas7bdat") library(haven) binary <- read_sas('binary.sas7bdat') names(binary) <- tolower(names(binary)) glm(admit~gre+gpa+as.factor(rank), family=binomial, data=binary) ## ## Call: glm(formula = admit ~ gre + gpa + as.factor(rank), family = binomial, ## data = binary) ## ## Coefficients: ## (Intercept) gre gpa ## -3.989979 0.002264 0.804038 ## as.factor(rank)2 as.factor(rank)3 as.factor(rank)4 ## -0.675443 -1.340204 -1.551464 ## ## Degrees of Freedom: 399 Total (i.e. Null); 394 Residual ## Null Deviance: 500 ## Residual Deviance: 458.5 AIC: 470.5 mle.logreg = function(fmla, data) { # Define the negative log likelihood function logl <- function(theta,x,y){ y <- y x <- as.matrix(x) beta <- theta[1:ncol(x)] # Use the log-likelihood of the Bernouilli distribution, where p is # defined as the logistic transformation of a linear combination # of predictors, according to logit(p)=(x%*%beta) loglik <- sum(-y*log(1 + exp(-(x%*%beta))) - (1-y)*log(1 + exp(x%*%beta))) return(-loglik) } # Prepare the data outcome = rownames(attr(terms(fmla),"factors"))[1] dfrTmp = model.frame(data) x = as.matrix(model.matrix(fmla, data=dfrTmp)) y = as.numeric(as.matrix(data[,match(outcome,colnames(data))])) # Define initial values for the parameters theta.start = rep(0,(dim(x)[2])) names(theta.start) = colnames(x) # Calculate the maximum likelihood mle = optim(theta.start,logl,x=x,y=y,hessian=F) # Obtain regression coefficients beta = mle$par # Calculate the Information matrix # The variance of a Bernouilli distribution is given by p(1-p) p = 1/(1+exp(-x%*%beta)) V = array(0,dim=c(dim(x)[1],dim(x)[1])) diag(V) = p*(1-p) IB = t(x)%*%V%*%x # Return estimates out = list(beta=beta,vcov=solve(IB),dev=2*mle$value) } binary$rank = factor(binary$rank) #Treat rank as a categorical variable fmla = as.formula("admit~gre+gpa+rank") #Create model formula mylogit = mle.logreg(fmla, binary) #Estimate coefficients mylogit ## $beta ## (Intercept) gre gpa rank2 rank3 ## -1.452667967 0.001346342 -0.130762174 0.719615014 -0.217101795 ## rank4 ## -0.010104814 ## ## $vcov ## (Intercept) gre gpa rank2 ## (Intercept) 1.2042376637 -2.685639e-04 -0.2751006299 -1.187537e-01 ## gre -0.0002685639 1.139736e-06 -0.0001266842 1.221842e-05 ## gpa -0.2751006299 -1.266842e-04 0.1022979333 7.058475e-03 ## rank2 -0.1187537263 1.221842e-05 0.0070584755 1.152379e-01 ## rank3 -0.1007985536 3.832995e-05 -0.0028178585 8.732829e-02 ## rank4 -0.1362108989 2.956123e-05 0.0089986268 8.826868e-02 ## rank3 rank4 ## (Intercept) -1.007986e-01 -1.362109e-01 ## gre 3.832995e-05 2.956123e-05 ## gpa -2.817858e-03 8.998627e-03 ## rank2 8.732829e-02 8.826868e-02 ## rank3 1.384753e-01 8.797479e-02 ## rank4 8.797479e-02 1.701939e-01 ## ## $dev ## [1] 497.3367 Previous Next