Up to 1.22

map handles linear models now

DESCRIPTION

@@ -1,19 +1,12 @@
 Package: rethinking
 Type: Package
 Title: Statistical Rethinking book package
-Version: 1.11
-Date: 2012-09-27
+Version: 1.22
+Date: 2013-09-26
 Author: Richard McElreath
 Maintainer: Richard McElreath <[email protected]>
 Depends:
-    MASS,
-    bbmle,
-    coda
+    MASS
+Imports: coda
 Description: Utilities for fitting and comparing models
-<<<<<<< HEAD
-License: GPLv3
-=======
-License: GPL version 2 or newer
-Collate:
-    'rethinking_utils.r'
->>>>>>> 2736e5bdc909542a79e08ad6922b529781c864fd
+License: GPLv3

NAMESPACE

@@ -0,0 +1,3 @@
+exportPattern(".")
+importFrom(coda, HPDinterval)
+importFrom(coda, as.mcmc)

R/coeftab.r

@@ -0,0 +1,108 @@
+# coeftab
+
+# coeftab class definition and show method
+setClass( "coeftab" , representation( coefs="matrix" , nobs="numeric" , AIC="numeric" , digits="numeric" , width="numeric" ) )
+coeftab.show <- function( object ) {
+    result <- object@coefs
+    if ( !is.null(object@nobs) ) {
+        result <- rbind( result , object@nobs )
+        rownames(result)[ nrow(result) ] <- "nobs"
+    }
+    coefs <- rrformat( result , digits=object@digits , width=object@width )
+    print( coefs , quote=FALSE , justify="right" )
+}
+setMethod( "show" , "coeftab" , function(object) coeftab.show(object) )
+
+coeftab <- function( ... , se=FALSE , se.inside=FALSE , nobs=TRUE , digits=2 , width=7 , rotate=FALSE , compare=FALSE ) {
+
+    # se=TRUE outputs standard errors
+    # se.inside=TRUE prints standard errors in parentheses in same column as estimates
+
+    if ( se.inside==TRUE ) se <- TRUE
+
+    # retrieve list of models
+    L <- list(...)
+    if ( is.list(L[[1]]) && length(L)==1 )
+        L <- L[[1]]
+
+    # retrieve model names from function call
+    mnames <- match.call()
+    mnames <- as.character(mnames)[2:(length(L)+1)]
+
+    # count number of unique parameters
+    param.names <- {}
+    for ( i in 1:length(L) ) {
+        c.names <- names( xcoef( L[[i]] ) )
+        param.names <- unique( c( param.names , c.names ) )
+    }
+    # columns for standard errors
+    if ( se==TRUE && se.inside==FALSE ) {
+        for ( i in 1:length(L) ) {
+            kse.names <- paste( names( xcoef( L[[i]] ) ) , ".se" , sep="" )
+            param.names <- unique( c( param.names , kse.names ) )
+        }
+    }
+
+    # make empty table
+    nk <- length(param.names)
+    d <- matrix( NA , ncol=nk )
+    d <- data.frame(d)
+    colnames(d) <- c( param.names )
+
+    # loop over models and insert values
+    for ( i in 1:length(L) ) {
+        klist <- xcoef( L[[i]] )
+        for ( j in 1:length(klist) ) {
+            d[i,][ names( klist[j] ) ] <- as.numeric( round( klist[j] , digits ) )
+        }
+    }
+    # insert standard errors
+    if ( se==TRUE ) {
+        for ( i in 1:length(L) ) {
+            kse <- xse( L[[i]] )
+            names(kse) <- names( xcoef( L[[i]] ) )
+            for ( j in 1:length(kse) ) {
+                if ( se.inside==FALSE )
+                    # own column
+                    d[i,][ paste(names(kse)[j],".se",sep="") ] <- as.numeric( round(kse[j],digits) )
+                else
+                    # combine with estimate
+                    d[i,][ names(kse)[j] ] <- paste( formatC( as.real(d[i,][ names(kse)[j] ]) , digits=digits ) , " (" , formatC( as.real( kse[j] ) , digits=digits ) , ")" , sep="" )
+            }
+        }
+    }
+
+    # add model names to rows
+    rownames(d) <- mnames
+
+    # formatting for parenthetical standard errors
+    if ( se.inside==TRUE && se==TRUE ) {
+        comment(d) <- "Values in parentheses are quadratic estimate standard errors."
+        colnames(d) <- paste( colnames(d) , "(se)" )
+        for ( i in 1:nrow(d) ) {
+            for ( j in 1:ncol(d) ) {
+                d[i,j] <- ifelse( is.na(d[i,j]) , "" , d[i,j] )
+            }
+        }
+    }
+
+    # add nobs
+    if ( nobs==TRUE ) {
+        nobs <- sapply( L , xnobs )
+    } else {
+        nobs <- 0
+    }
+
+    # add AICc and weights
+    if ( compare==TRUE ) {
+        require(bbmle)
+        d$AICc <- sapply( 1:length(L) , function(z) AICc( L[[z]] , nobs=d$nobs[z] ) )
+        deltAICc <- d$AICc - min(d$AICc) 
+        wAIC <- exp( -0.5*deltAICc )
+        d$weight <- wAIC / sum( wAIC )
+    }
+
+    # return table
+    if ( rotate==FALSE ) d <- t(d) # models along top is default
+    new( "coeftab" , coefs=d , nobs=nobs , digits=digits , width=width )
+}

R/colors.r

@@ -0,0 +1,20 @@
+# color utility functions
+
+col.desat <- function( acol , amt=0.5 ) {
+    acol <- col2rgb(acol)
+    ahsv <- rgb2hsv(acol)
+    ahsv[2] <- ahsv[2] * amt
+    hsv( ahsv[1] , ahsv[2] , ahsv[3] )
+}
+
+col.alpha <- function( acol , alpha=0.2 ) {
+    acol <- col2rgb(acol)
+    acol <- rgb(acol[1]/255,acol[2]/255,acol[3]/255,alpha)
+    acol
+}
+
+col.dist <- function( x , mu=0 , sd=1 , col="slateblue" ) {
+    cols <- sapply( x , function(z) exp(-(z-mu)^2/sd) )
+    cols <- sapply( cols , function(z) col.alpha(col,z) )
+    cols
+}

R/compare.r

@@ -0,0 +1,66 @@
+# compare
+
+# compare class definition and show method
+setClass( "compareIC" , representation( output="data.frame" ) )
+compare.show <- function( object ) {
+    print( round( object@output , 2 ) )
+}
+setMethod( "show" , "compareIC" , function(object) compare.show(object) )
+
+# AICc/BIC model comparison table
+compare <- function( ... , nobs=NULL , sort="AICc" , delta=TRUE , digits=4 ) {
+    require(bbmle)
+
+    if ( is.null(nobs) ) {
+        stop( "Must specify number of observations (nobs)." )
+    }
+
+    getdf <- function(x) {
+        if (!is.null(df <- attr(x, "df"))) 
+            return(df)
+        else if (!is.null(df <- attr(logLik(x), "df"))) 
+            return(df)
+    }
+
+    # need own BIC, as one in stats doesn't allow nobs
+    myBIC <- function(x,nobs) {
+        k <- getdf(x)
+        as.numeric( -2*logLik(x) + log(nobs)*k )
+    }
+
+    # retrieve list of models
+    L <- list(...)
+    if ( is.list(L[[1]]) && length(L)==1 )
+        L <- L[[1]]
+
+    # retrieve model names from function call
+    mnames <- match.call()
+    mnames <- as.character(mnames)[2:(length(L)+1)]
+
+    AICc.list <- sapply( L , function(z) AICc( z , nobs=nobs ) )
+    BIC.list <- sapply( L , function(z) myBIC( z , nobs=nobs ) )
+    dAICc <- AICc.list - min( AICc.list )
+    dBIC <- BIC.list - min( BIC.list )
+
+    post.AICc <- exp( -0.5*dAICc ) / sum( exp(-0.5*dAICc) )
+    post.BIC <- exp( -0.5*dBIC ) / sum( exp(-0.5*dBIC) )
+
+    #post.AICc <- format.pval( post.AICc , digits=digits )
+    #post.BIC <- format.pval( post.BIC , digits=digits )
+
+    k <- sapply( L , getdf )
+
+    result <- data.frame( k=k , AICc=AICc.list , BIC=BIC.list , w.AICc=post.AICc , w.BIC=post.BIC )
+    if ( delta==TRUE ) {
+        r2 <- data.frame( dAICc=dAICc , dBIC=dBIC )
+        result <- cbind( result , r2 )
+    }
+    rownames( result ) <- mnames
+
+    if ( !is.null(sort) ) {
+        result <- result[ order( result$AICc ) , ]
+    }
+
+    new( "compareIC" , output=result )
+}
+

R/distributions.r

@@ -0,0 +1,138 @@
+# distributions
+
+# convert to sigma from tau
+tau2sig <- function( tau ) 1/abs(tau)
+
+# Gaussian density that takes tau instead of sigma
+dnormtau <- function( x , mean=0 , tau=1 , log=FALSE ) dnorm( x , mean=mean , sd=tau2sig( tau ) , log=log )
+
+rnormtau <- function( n , mean=0 , tau=1 ) rnorm( n=n , mean=mean , sd=tau2sig(tau) )
+
+# ordered categorical density functions
+logistic <- function( x ) {
+    p <- 1 / ( 1 + exp( -x ) )
+    p <- ifelse( x==Inf , 1 , p )
+    p
+}
+
+pordlogit <- function( x , a , phi , log=FALSE ) {
+    a <- c( as.numeric(a) , Inf )
+    if ( length(phi) == 1 ) {
+        p <- logistic( a[x] - phi )
+    } else {
+        p <- matrix( NA , ncol=length(x) , nrow=length(phi) )
+        for ( i in 1:length(phi) ) {
+            p[i,] <- logistic( a[x] - phi[i] )
+        }
+    }
+    if ( log==TRUE ) p <- log(p)
+    p
+}
+
+dordlogit <- function( x , a , phi , log=FALSE ) {
+    a <- c( as.numeric(a) , Inf )
+    p <- logistic( a[x] - phi )
+    na <- c( -Inf , a )
+    np <- logistic( na[x] - phi )
+    p <- p - np
+    if ( log==TRUE ) p <- log(p)
+    p
+}
+
+rordlogit <- function( n , a , phi=0 ) {
+    a <- c( as.numeric(a) , Inf )
+    k <- 1:length(a)
+    p <- dordlogit( k , a=a , phi=phi , log=FALSE )
+    y <- sample( k , size=n , replace=TRUE , prob=p )
+    y
+}
+
+# beta density functions parameterized as prob,theta
+
+dbeta2 <- function( x , prob , theta , log=FALSE ) {
+    a <- prob * theta
+    b <- (1-prob) * theta
+    dbeta( x , shape1=a , shape2=b , log=log )
+}
+
+rbeta2 <- function( n , prob , theta ) {
+    a <- prob * theta
+    b <- (1-prob) * theta
+    rbeta( n , shape1=a , shape2=b )
+}
+
+# based on Ben Bolker's dbetabinom in emdbook package
+dbetabinom <- function (x, prob, size, theta, shape1, shape2, log = FALSE) 
+{
+    if (missing(prob) && !missing(shape1) && !missing(shape2)) {
+        prob <- shape1/(shape1 + shape2)
+        theta <- shape1 + shape2
+    }
+    v <- lfactorial(size) - lfactorial(x) - lfactorial(size - 
+        x) - lbeta(theta * (1 - prob), theta * prob) + lbeta(size - 
+        x + theta * (1 - prob), x + theta * prob)
+    if (any(n <- nonint(x))) {
+        warning("non-integer x detected; returning zero probability")
+        v[n] <- -Inf
+    }
+    if (log) 
+        v
+    else exp(v)
+}
+
+# gamma-poisson functions
+
+dgamma2 <- function( x , mu , scale , log=FALSE ) {
+    dgamma( x , shape=mu / scale , scale=scale , log=log )
+}
+
+rgamma2 <- function( n , mu , scale ) {
+    rgamma( n , shape=mu/scale , scale=scale )
+}
+
+dgampois <- function( x , mu=NULL , shape=NULL , scale=NULL , rate=NULL , log=FALSE ) {
+    if ( !is.null(rate) ) scale <- 1/rate
+    if ( !is.null(mu) ) shape <- mu / scale
+    prob <- 1 / ( 1 + scale )
+    dnbinom( x , size=shape , prob=prob , log=log )
+}
+
+rgampois <- function( n , mu=NULL , shape=NULL , scale=NULL , rate=NULL ) {
+    if ( !is.null(rate) ) scale <- 1/rate
+    if ( !is.null(mu) ) shape <- mu / scale
+    prob <- 1 / ( 1 + scale )
+    rnbinom( n , size=shape , prob=prob )
+}
+
+###
+# zero-inflated gamma functions
+dzerogamma <- function( x , prob , mu , scale , log=TRUE ) {
+    p1 <- dgamma2( x , mu=mu , s=scale , log=FALSE )
+    p <- ifelse( x==0 , prob , (1-prob)*p1 )
+    if ( log==TRUE ) p <- log(p)
+    p
+}
+
+# laplace (double exponential)
+dlaplace <- function(x,lambda,location=0,log=FALSE) {
+    # f(y) = (1/(2b)) exp( -|y-a|/b )
+    l <- (1/(2*lambda))*exp( -abs(x-location)/lambda )
+    if ( log==TRUE ) l <- log(l)
+    l
+}
+
+rlaplace <- function(n,lambda,location=0) {
+    # generate an exponential sample, then generate a random sign, then add location
+    y <- rexp(n=n,rate=lambda) * sample( c(-1,1) , size=n , replace=TRUE )
+    y <- y + location
+    y
+}
+
+### softmax rule for multinomial
+multilogistic <- function( x , lambda=1 , diff=TRUE , log=FALSE ) {
+    if ( diff==TRUE ) x <- x - min(x)
+    f <- exp( lambda * x )
+    if ( log==FALSE ) result <- f / sum(f)
+    if ( log==TRUE ) result <- log(f) - log(sum(f))
+    result
+}