R version 2.15.1 (2012-06-22) -- "Roasted Marshmallows" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: x86_64-pc-mingw32/x64 (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > pkgname <- "LogReg" > source(file.path(R.home("share"), "R", "examples-header.R")) > options(warn = 1) > options(pager = "console") > library('LogReg') Loading required package: MASS > > assign(".oldSearch", search(), pos = 'CheckExEnv') > cleanEx() > nameEx("est.bes") > ### * est.bes > > flush(stderr()); flush(stdout()) > > ### Name: est.bes > ### Title: A function to estimate betas in a logistic regression > ### Aliases: est.bes > ### Keywords: ~kwd1 ~kwd2 > > ### ** Examples > > ##---- Should be DIRECTLY executable !! ---- > ##-- ==> Define data, use random, > ##-- or do help(data=index) for the standard data sets. > > ## The function is currently defined as > function (Y, X, intercept = TRUE, start.vals = rep(0, ncol(X) + + as.numeric(intercept))) + { + x <- start.vals + opt <- optim(x, fn = lr.loglike, gr = NULL, Y = Y, X = X, + intercept = intercept, method = "BFGS", control = list(fnscale = -1), + hessian = T) + beta.hat <- opt$par + beta.cov <- solve(opt$hessian)/length(Y) + if (intercept) { + Xm <- cbind(rep(1, length(Y)), X) + } + else { + Xm <- X + } + fitted.lp <- Xm %*% beta.hat + fitted.Y <- expit(fitted.lp) + conv <- TRUE + if (opt$convergence != 0) { + conv <- FALSE + print("Optim did not converge") + } + return(list(ests = beta.hat, varcov = beta.cov, fitted.lp = fitted.lp, + fitted.Y = fitted.Y, conv = conv)) + } function (Y, X, intercept = TRUE, start.vals = rep(0, ncol(X) + as.numeric(intercept))) { x <- start.vals opt <- optim(x, fn = lr.loglike, gr = NULL, Y = Y, X = X, intercept = intercept, method = "BFGS", control = list(fnscale = -1), hessian = T) beta.hat <- opt$par beta.cov <- solve(opt$hessian)/length(Y) if (intercept) { Xm <- cbind(rep(1, length(Y)), X) } else { Xm <- X } fitted.lp <- Xm %*% beta.hat fitted.Y <- expit(fitted.lp) conv <- TRUE if (opt$convergence != 0) { conv <- FALSE print("Optim did not converge") } return(list(ests = beta.hat, varcov = beta.cov, fitted.lp = fitted.lp, fitted.Y = fitted.Y, conv = conv)) } > > > > cleanEx() > nameEx("expit") > ### * expit > > flush(stderr()); flush(stdout()) > > ### Name: expit > ### Title: The Expit Function > ### Aliases: expit > ### Keywords: ~kwd1 ~kwd2 > > ### ** Examples > > ##---- Should be DIRECTLY executable !! ---- > ##-- ==> Define data, use random, > ##-- or do help(data=index) for the standard data sets. > > ## The function is currently defined as > function (x) + { + 1/(1 + exp(-x)) + } function (x) { 1/(1 + exp(-x)) } > > > > cleanEx() > nameEx("gbd") > ### * gbd > > flush(stderr()); flush(stdout()) > > ### Name: gbd > ### Title: A function to generate some binary regression data > ### Aliases: gbd > ### Keywords: ~kwd1 ~kwd2 > > ### ** Examples > > ##---- Should be DIRECTLY executable !! ---- > ##-- ==> Define data, use random, > ##-- or do help(data=index) for the standard data sets. > > ## The function is currently defined as > function (n = 1000, p = 1, intercept = TRUE, betas = rep(0, p + + intercept), mu = rep(0, p), Sigma = diag(1, p), seed = 1) + { + require(MASS) + set.seed(seed) + X <- mvrnorm(n, mu, Sigma) + Xm <- X + if (intercept) { + Xm <- cbind(rep(1, n), X) + } + lp <- Xm %*% betas + pv <- expit(lp) + Y <- rbinom(n, 1, pv) + return(list(Y = Y, X = X)) + } function (n = 1000, p = 1, intercept = TRUE, betas = rep(0, p + intercept), mu = rep(0, p), Sigma = diag(1, p), seed = 1) { require(MASS) set.seed(seed) X <- mvrnorm(n, mu, Sigma) Xm <- X if (intercept) { Xm <- cbind(rep(1, n), X) } lp <- Xm %*% betas pv <- expit(lp) Y <- rbinom(n, 1, pv) return(list(Y = Y, X = X)) } > > > > cleanEx() > nameEx("logit") > ### * logit > > flush(stderr()); flush(stdout()) > > ### Name: logit > ### Title: The logit function > ### Aliases: logit > ### Keywords: ~kwd1 ~kwd2 > > ### ** Examples > > ##---- Should be DIRECTLY executable !! ---- > ##-- ==> Define data, use random, > ##-- or do help(data=index) for the standard data sets. > > ## The function is currently defined as > function (x) + { + return(log(x/(1 - x))) + } function (x) { return(log(x/(1 - x))) } > > > > cleanEx() > nameEx("lr.loglike") > ### * lr.loglike > > flush(stderr()); flush(stdout()) > > ### Name: lr.loglike > ### Title: Evaluate the logistic regression log-likelihood > ### Aliases: lr.loglike > ### Keywords: ~kwd1 ~kwd2 > > ### ** Examples > > ##---- Should be DIRECTLY executable !! ---- > ##-- ==> Define data, use random, > ##-- or do help(data=index) for the standard data sets. > > ## The function is currently defined as > function (betas, Y, X, intercept = TRUE) + { + if (!intercept) { + lp <- X %*% betas + } + else { + lp <- cbind(rep(1, length(Y)), X) %*% betas + } + pv <- expit(lp) + return(sum(dbinom(Y, 1, pv, log = T))) + } function (betas, Y, X, intercept = TRUE) { if (!intercept) { lp <- X %*% betas } else { lp <- cbind(rep(1, length(Y)), X) %*% betas } pv <- expit(lp) return(sum(dbinom(Y, 1, pv, log = T))) } > > > > cleanEx() > nameEx("plot.lrsurf") > ### * plot.lrsurf > > flush(stderr()); flush(stdout()) > > ### Name: plot.lrsurf > ### Title: Plot the log-likelihood surface for a single predictor > ### Aliases: plot.lrsurf > ### Keywords: ~kwd1 ~kwd2 > > ### ** Examples > > ##---- Should be DIRECTLY executable !! ---- > ##-- ==> Define data, use random, > ##-- or do help(data=index) for the standard data sets. > > ## The function is currently defined as > function (Y, X, b0 = seq(-2, 2, by = 0.01), b1 = seq(-2, 2, by = 0.01), + ncl = 20) + { + lrs <- outer(b0, b1, "+") * 0 + for (i in 1:length(Y)) { + lpv <- outer(b0, b1 * X[i], "+") + lrs <- lrs + lpv * Y[i] - log(1 + exp(lpv)) + } + image(b0, b1, lrs, xlab = expression(beta[0]), ylab = expression(beta[1])) + contour(b0, b1, lrs, add = T, nlevels = ncl) + } function (Y, X, b0 = seq(-2, 2, by = 0.01), b1 = seq(-2, 2, by = 0.01), ncl = 20) { lrs <- outer(b0, b1, "+") * 0 for (i in 1:length(Y)) { lpv <- outer(b0, b1 * X[i], "+") lrs <- lrs + lpv * Y[i] - log(1 + exp(lpv)) } image(b0, b1, lrs, xlab = expression(beta[0]), ylab = expression(beta[1])) contour(b0, b1, lrs, add = T, nlevels = ncl) } > > > > cleanEx() > nameEx("test.data") > ### * test.data > > flush(stderr()); flush(stdout()) > > ### Name: test.data > ### Title: A test binary data set with two predictors > ### Aliases: test.data > ### Keywords: datasets > > ### ** Examples > > data(test.data) > ## maybe str(test.data) ; plot(test.data) ... > > > > ### *