|
|
ordinarylogit <- function(b, X, Y){
X <- cbind(1, X) #add constant
xb <- X %*% b
llik <- (-Y * log(1 + exp(-xb))) - (1 - Y) * log(1 + exp(xb)) #the actual log-likelihood
sum(llik) #pass results back to optim
}
|
- Cloglog likelihood routine
|
cloglog <- function(b, X, Y) {
xb <- X %*% b
llik <- Y * log(1 - exp(-exp(xb))) - (1 - Y) * (exp(xb))
sum(llik)
}
|
- Typical maximization call using the optim command
|
stval <- runif(ncol(indepvar.com) + 1)*0.01 #randomly generate starting values for all covariates plus constant
mymodel <- optim(stval, cloglog, hessian = TRUE, method = "BFGS", control = list(fnscale = -1, reltol = 1e-20, maxit = 1000), X = indepvar.com, Y = depvar)
|
- Table that summarizes estimation results for models without constant
|
optim.table.par0 <- function(results, varnames = NULL, tit = NULL) {
parameter <- as.matrix(results$par)
vacov <- as.matrix(solve(-1 * results$hessian))
stderror <- as.matrix(sqrt(diag(vacov)))
tstat <- as.matrix(parameter / stderror)
tsub1 <- cbind(round(parameter, digits = 3), round(stderror, digits = 3), round(tstat, digits = 3))
table1 <- cbind(tsub1, round(2 * pnorm(abs(tstat), lower.tail = FALSE), digits=3))
dimnames(table1) <- list(varnames, c("Estimate", "Std. Error", "z-value", "p-value"))
cat("\n\n\n ", tit, "\n\n")
print(table1)
}
|
- Table that summarizes estimation results for models with constant
|
optim.table.par1 <- function(results, varnames=NULL, tit=NULL) {
varnames <- c("Constant", varnames[1:length(varnames)])
parameter <- as.matrix(results$par)
vacov <- as.matrix(solve(-1 * results$hessian))
stderror <- as.matrix(sqrt(diag(vacov)))
tstat <- as.matrix(parameter / stderror)
tsub1 <- cbind(round(parameter, digits=3), round(stderror, digits=3), round(tstat, digits=3))
table1 <- cbind(tsub1, round(2*pnorm(abs(tstat), lower.tail=FALSE), digits=3))
dimnames(table1)<-list(varnames, c("Estimate", "Std. Error", "z-value", "p-value"))
print(table1)
}
|
- The likelihood-ratio test
|
lr.test <- function(restricted, full){
dgf <- abs(length(restricted$par)-length(full$par))
r <- -2 * (restricted$value - full$value)
p <- 1-pchisq(r, dgf)
cat("degrees of freedom: ", dgf, "\n")
cat("likelihood ratio statistic: ", r,"\n")
cat("p value: ", p, " (under H0 that the restricted model is correct)\n\n")
list(r = r, p = p)
}
|
|
|
vuong <- function(y, x1, x2, mod1, mod2, level = 0.05){
LR <- mod1$value - mod2$value #likelihood ratio
xb1 <- as.matrix(x1) %*% as.matrix(mod1$par)
xb2 <- as.matrix(x2) %*% as.matrix(mod2$par)
llik <- y*log(1-exp(-exp(xb1)))-(1-y)*(exp(xb1)) -
y*log(1-exp(-exp(xb2)))-(1-y)*(exp(xb2))
sum.llik <- sum(llik)
llik.sq <- llik^2
omega.sq <- 1 / length(y) * sum(llik.sq) - (1 / length(y) * sum.llik)^2
V <- 1 / sqrt(omega.sq * length(y)) * LR
if (V > qnorm(1-(level/2)))cat("Model 1 is superior \n\n")
if (V < qnorm(level/2)) cat("Model 2 is superior \n\n")
if (abs(V) <= qnorm(1-(level/2))) cat("The models are not distinguishable \n\n")
p <- pnorm(V)
list(vstat = V, pvalue = 1-p)
}
|
- Plotting with control over lines types, labels, and legend.
|
plot(c(1,20),c(0,1),main=" Title",xlab="",ylab="", type="n")
lines(obj1, lty=1)
lines(obj2, lty=3)
lines(obj3, lty=2)
lines(obj4, lty=4)
legend(1,1, legend=c("cat. 1","cat. 2", "cat. 3", "cat. 4"), lty=c(1,2,3,4 ))
|
|
|
|
|
|
|