##------------------------------------------------------------##
##  Script for Object-Oriented R Programming & Lexical Scope  ##
##                        John Fox                            ##
##                    Programming in R                        ##
##                ICPSR Summer Program 2009                   ##
##------------------------------------------------------------##

# The S3 object system

    # S3 classes

library(car) # for data

mod.prestige <- lm(prestige ~ income + education + women,
    data=Prestige)
attributes(mod.prestige)
class(mod.prestige)


v <- 1:10
v
attributes(v)
class(v)

class(v) <- "character"
attributes(v)
v
class(v)

    # S3 generic functions and methods
    
print # the print generic

print.lm # print method for "lm" objects

mod.prestige
print(mod.prestige) # equivalent
print.lm(mod.prestige)  # equivalent, but bad form

methods("print") # print methods
methods(class="lm") # methods for objects of class "lm"

    # S3 "inheritance"

mod.mroz <- glm(lfp ~ ., family=binomial, data=Mroz)
class(mod.mroz)

    # Example: a logistic-regression function

lreg3 <- function(X, y, predictors=colnames(X), max.iter=10,
        tol=1E-6, constant=TRUE)
{
    if (!is.numeric(X) || !is.matrix(X))  # data checks
        stop("X must be a numeric matrix")
    if (!is.numeric(y) || !all(y == 0 | y == 1))
        stop("y must contain only 0s and 1s")
    if (nrow(X) != length(y))
        stop("X and y contain different numbers of observations")
    if (constant) { # attach constant?
        X <- cbind(1, X)
        colnames(X)[1] <- "Constant"
    }
    b <- b.last <- rep(0, ncol(X))
    it <- 1
    while (it <= max.iter){
        p <- as.vector(1/(1 + exp(-X %*% b)))
        var.b <- solve(crossprod(X, p * (1 - p) * X))
        b <- b + var.b %*% crossprod(X, y - p)
        if (max(abs(b - b.last)/(abs(b.last) + 0.01*tol)) < tol) break
        b.last <- b
        it <- it + 1
    }
    if (it > max.iter) warning("maximum iterations exceeded")
    dev <- -2*sum(y*log(p) + (1 - y)*log(1 - p))
    result <- list(coefficients=as.vector(b), var=var.b,
        deviance=dev, converged= it <= max.iter,
        predictors=predictors)
    class(result) <- "lreg3"  # assign class
    result
}

Mroz$lfp <- with(Mroz, ifelse(lfp == "yes", 1, 0))
Mroz$wc <- with(Mroz, ifelse(wc == "yes", 1, 0))
Mroz$hc <- with(Mroz, ifelse(hc == "yes", 1, 0))

mod.mroz.3 <- with(Mroz, lreg3(cbind(k5, k618, age, wc, hc, lwg, inc), lfp))
class(mod.mroz.3)
mod.mroz.3 # whoops!

print.lreg3 <- function(x, ...)  # print method for class "lreg3"
{
    coef <- x$coefficients
    names(coef) <- x$predictors
    print(coef)
    if (!x$converged) cat("\n *** lreg did not converge ***\n")
    invisible(x) # note: passes through argument invisibly
}

mod.mroz.3

summary  # summary generic

summary.lreg3 <- function(object, ...)  # summary method for class "lreg3"
{
    b <- object$coefficients
    se <- sqrt(diag(object$var))
    z <- b/se
    table <- cbind(b, se, z, 2*(1-pnorm(abs(z))))
    colnames(table) <- c("Estimate", "Std.Err", "Z value", "Pr(>z)")
    rownames(table) <- object$predictors
    result <- list(coef=table, deviance=object$deviance,
        converged=object$converged)
    class(result) <- "summary.lreg3" # creates an object of class "summary.lreg3"
    result
}

print.summary.lreg3 <- function(x, ...) # print method for class "summary.lreg3"
{
    printCoefmat(x$coef)
    cat("\nDeviance =", x$deviance,"\n")
    if (!x$converged) cat("\n Note: *** lreg did not converge ***\n")
}

summary(mod.mroz.3)

    # writing a generic function
    
names(summary(mod.prestige))
    
rsq <- function(model, ...)
{
    UseMethod("rsq")
}

rsq.lm <- function(model, adjusted=FALSE, ...)
{
    summary <- summary(model)
    if (adjusted) summary$adj.r.squared else summary$r.squared
}

rsq(mod.prestige)
rsq(mod.prestige, adjusted=TRUE)

rsq(mod.mroz) #  via inheritance (doesn't work)

names(summary(mod.mroz))
names(mod.mroz)

rsq.glm <- function(model, ...) 
{
    1 - model$deviance/model$null.deviance
}

rsq(mod.mroz)

# The S4 Object System

    # definition of S4 classes
    
setClass("lreg4",
    representation(coefficients="numeric", var="matrix",
        iterations="numeric", deviance="numeric",
        predictors="character"))
        
lreg4 <- function(X, y, predictors=colnames(X), constant=TRUE,
max.iter=10, tol=1E-6)
{
    if (!is.numeric(X) || !is.matrix(X))
        stop("X must be a numeric matrix")
    if (!is.numeric(y) || !all(y == 0 | y == 1))
        stop("y must contain only 0s and 1s")
    if (nrow(X) != length(y))
        stop("X and y contain different numbers of observations")
    if (constant) {
        X <- cbind(1, X)
        colnames(X)[1] <- "Constant"
    }
    b <- b.last <- rep(0, ncol(X))
    it <- 1
    while (it <= max.iter){
        p <- as.vector(1/(1 + exp(-X %*% b)))
        var.b <- solve(crossprod(X, p * (1 - p) * X))
        b <- b + var.b %*% crossprod(X, y - p)
        if (max(abs(b - b.last)/(abs(b.last) + 0.01*tol)) < tol) break
        b.last <- b
        it <- it + 1
    }
    if (it > max.iter) warning("maximum iterations exceeded")
    # create an instance of the "lreg4" class:
    result <- new("lreg4", coefficients=as.vector(b),
        var=var.b, iterations=it,
        deviance=-2*sum(y*log(p) + (1 - y)*log(1 - p)),
        predictors=predictors)
    result
}

mod.mroz.4 <- with(Mroz, lreg4(cbind(k5, k618, age, wc, hc, lwg, inc), lfp))
class(mod.mroz.4)
mod.mroz.4

show  # the S4 generic function show

    # defining an S4 method
    
setMethod("show", signature(object="lreg4"),
    definition=function(object)
        {
            coef <- object@coefficients
            names(coef) <- object@predictors
            print(coef)
        }
    )

mod.mroz.4  # invokes show method

setMethod("summary", signature(object="lreg4"),
    definition=function(object, ...)
        {
            b <- object@coefficients
            se <- sqrt(diag(object@var))
            z <- b/se
            table <- cbind(b, se, z, 2*(1-pnorm(abs(z))))
            colnames(table) <- c("Estimate", "Std.Err", "Z value", "Pr(>z)")
            rownames(table) <- object@predictors
            printCoefmat(table)
            cat("\nDeviance =", object@deviance,"\n")
        }
    )
    
summary(mod.mroz.4)


# Lexical scope

f <- function (x) x + a
a <- 10
x <- 5
f(2) # x bound to 2 in frame of f(), a to 10 in global frame
x    # global x is undisturbed

f <- function (x)
{
    a <- 5
    g(x)
}
g <- function(y) y + a
f(2) # a bound to 10 in global frame
a    # global a is undisturbed


f <- function (x)
{
    a <- 5
    g <- function (y) y + a
    g(x)
}
f(2) # a is bound to 5, x to 2 in frame of f(), y to 2 in frame of g()

    # a function that returns a closure (function + environment)
    
makePower <- function(power)
{
    function(x) x^power
}

square <- makePower(2)
square # power bound to 2
square(4)

cuberoot <- makePower(1/3)
cuberoot # power bound to 1/3
cuberoot(64)