# negative binomial data
prep.bar <- function( d , dmax=max(d) ) {
    x <- rep(0,dmax+1)
    for ( i in 1:(dmax+1) ) {
        x[i] <- length( d[d==(i-1)] )
    }
    x
}

par(mfcol=c(2,2))

y <- rnbinom( 1000 , mu=2 , size=10 )
barplot( prep.bar(y) , names.arg=0:max(y) , xlab="y" , ylab="Frequency" , main="mu=2,size=10" ,xlim=c(0,18) , space=0 )

y <- rnbinom( 1000 , mu=6 , size=10 )
barplot( prep.bar(y) , names.arg=0:max(y) , xlab="y" , ylab="Frequency" , main="mu=6,size=10" ,xlim=c(0,18) , space=0 )

y <- rnbinom( 1000 , mu=2 , size=2 )
barplot( prep.bar(y) , names.arg=0:max(y) , xlab="y" , ylab="Frequency" , main="mu=2,size=2" ,xlim=c(0,18) , space=0 )

y <- rnbinom( 1000 , mu=6 , size=2 )
barplot( prep.bar(y) , names.arg=0:max(y) , xlab="y" , ylab="Frequency" , main="mu=6,size=2" ,xlim=c(0,18) , space=0 )

# nbinom converges to poisson, as n goes to infinity
par(mfcol=c(2,3))

barplot( prep.bar(rnbinom(10000,mu=5,size=1))  , xlab="y" , ylab="Frequency" , space=0 , main="nbinom, size=1" )

barplot( prep.bar(rpois(10000,lambda=5))  , xlab="y" , ylab="Frequency" , space=0 , main="poisson, same mean")

barplot( prep.bar(rnbinom(10000,mu=5,size=10))  , xlab="y" , ylab="Frequency" , space=0 , main="nbinom, size=10")

barplot( prep.bar(rpois(10000,lambda=5))  , xlab="y" , ylab="Frequency" , space=0 , main="poisson, same mean")

barplot( prep.bar(rnbinom(10000,mu=5,size=1000)) , xlab="y" , ylab="Frequency" , space=0 , main="nbinom, size=1000")

barplot( prep.bar(rpois(10000,lambda=5))  , xlab="y" , ylab="Frequency" , space=0 , main="poisson, same mean")

# gamma mixing to get nbinom
my.dgamma <- function( x , k , t , log=FALSE ) {
    x^(k-1) * exp(-x/t) / ( t^k * gamma(k) )
}
curve( my.dgamma(x,3,2) , from=0 , to=20 , xlab="lambda (mean of poisson)" , ylab="probability" )
# equivalent to using built-in function:
# curve( dgamma(x,shape=3,scale=2) , from=0 , to=20 , xlab="lambda (mean of poisson)" , ylab="probability" )
lines( c(3*2,3*2) , c(0,my.dgamma(3*2,3,2)) , lty=2 )

y <- rpois( 100 , 1 )
barplot( prep.bar(y,dmax=20) , names.arg=0:20 , xlab="count" , ylab="Frequency" , xlim=c(0,20) , space=0 )

y <- rpois( 100 , 13 )
barplot( prep.bar(y,dmax=20) , names.arg=0:20 , xlab="count" , ylab="Frequency" , xlim=c(0,20) , space=0 )

y <- rpois( 100 , 4.5 )
barplot( prep.bar(y,dmax=20) , names.arg=0:20 , xlab="count" , ylab="Frequency" , xlim=c(0,20) , space=0 )

# demonstration of how gamma-poisson produces nbinom distribution
rgammapois <- function( n , k , t ) {
    lambdas <- rgamma( n , shape=k , scale=t )
    rpois( n , lambda=lambdas )
}
# y1 and y2 will converge
y1 <- rgammapois(100000,3,2)
y2 <- rnbinom(100000,size=3,mu=3*2) # size = k; mu = k*t
barplot( prep.bar(y1) , names.arg=0:max(y1) , xlab="count" , ylab="Frequency" , xlim=c(0,max(y1)) , space=0 )
barplot( prep.bar(y2) , names.arg=0:max(y2) , xlab="count" , ylab="Frequency" , xlim=c(0,max(y2)) , space=0 )

# interactive progressive sampling of gamma-poisson
n <- 200
k <- 3
th <- 2
par(mfcol=c(2,1),ask=FALSE)
curve( dgamma(x,shape=k,scale=th) , from=0 , to=20 , xlab="lambda (mean of poisson)" , ylab="probability" )
barplot( prep.bar(y2) , names.arg=0:max(y2) , xlab="count" , ylab="Frequency" , xlim=c(0,max(y2)) , space=0 , col="white" , border=NA )
y <- {}
l <- {}
for ( i in 1:n ) {
    rl <- rgamma( 1 , shape=k , scale=th )
    l <- c(l,rl)
    curve( dgamma(x,shape=k,scale=th) , from=0 , to=20 , xlab="lambda (mean of poisson)" , ylab="probability" , main="Gamma distribution" )
    for ( j in 1:length(l) ) lines( c(l[j],l[j]) , c(0, dgamma(l[j],shape=k,scale=th) ) )
    lines( c(rl,rl) , c(0, dgamma(rl,shape=k,scale=th) ) , col="red" , lwd=3 )
    y <- c( y , rpois(2,rl) )
    barplot( prep.bar(y,dmax=35) , names.arg=0:35 , xlab="count" , ylab="Frequency" , xlim=c(0,35) , space=0 , main="Poisson draws with Gamma means" )
    par(ask=TRUE)
}
par(ask=FALSE)

# visualize our sampled gamma-poisson against ideal poisson with same mean as gamma
y1 <- prep.bar(y,dmax=35)
barplot( y1 , names.arg=0:35 , xlab="count" , ylab="Frequency" , xlim=c(0,35) , space=0 , main="Poisson draws with Gamma means" )
y2 <- sapply( 0:35 , function(z) dpois(z,lambda=6) )
y2 <- max(y1) * y2 / max(y2)
lines( 0:35 , y2 , col="red" )

# back to salamanders
library(bbmle)
d <- read.csv("salamanders.csv")
barplot( prep.bar(d$SALAMAN) , names.arg=0:max(d$SALAMAN) , xlab="salamander density" , ylab="frequency" )

mp <- mle2( d$SALAMAN ~ dpois( lambda=a ) , start=list(a=mean(d$SALAMAN)) )
mnb <- mle2( d$SALAMAN ~ dnbinom( mu=a , size=exp(n) ) , start=list(a=mean(d$SALAMAN),n=0) )

# generate predicted frequencies for poisson model
predict.p <- dpois( 0:13 , coef(mp)[1] ) * length(d$SALAMAN)
# generate predicted frequencies for neg-binom model
predict.nb <- dnbinom( 0:13 , mu=coef(mnb)[1] , size=exp(coef(mnb)[2]) ) * length(d$SALAMAN)

# overlay predictions on data
barplot( prep.bar(d$SALAMAN) , names.arg=0:max(d$SALAMAN) , xlab="salamander density" , ylab="frequency" , space=0 )
lines( 0:13 , predict.p , col="red" )
lines( 0:13 , predict.nb , col="blue" )

# likelihood changes very little, as n approaches infinity. so need to start search from low values, or risk just hanging out in the flat space on the right of this plot...
x <- seq( from=0 , to=100 , by=0.1 )
liks <- sapply( x , function(z) -sum(dnbinom( d$SALAMAN , mu=coef(mnb)[1] , size=z , log=TRUE )) )
plot( x , liks , xlab="size (dispersion)" , ylab="-log(likelihood)" , type="l" )
# another view, even more extreme
x <- seq( from=0 , to=1000 , by=0.1 )
liks <- sapply( x , function(z) -sum(dnbinom( d$SALAMAN , mu=coef(mnb)[1] , size=z , log=TRUE )) )
plot( x , liks , xlab="size (dispersion)" , ylab="-log(likelihood)" , type="l" )

plot( x , liks , xlab="size (dispersion)" , ylab="-log(likelihood)" , type="l" , log="x" )

# now imagine fitting the parameters of the underlying gamma, instead
# then mu = k*th, size = k
mgp <- mle2( d$SALAMAN ~ dnbinom( mu=k*th , size=k ) , start=list(k=3,th=2) )
# estimates are equivalent to mnb fit
# plot underlying gamma, using either formulation
k <- coef(mgp)
curve( dgamma( x , shape=k[1] , scale=k[2] ) , from=0 , to=20 )
# looking at above, a very skewed distribution, where most sites have almost no salamanders in them. i.e. we are looking at the inferred heterogeneity among sites, when we look at the gamma dist above
# can also plot from straight mnb estimates:
k <- coef(mnb)
# now shape = k = exp(n) , scale = th = mu/th = a/exp(n)
curve( dgamma( x , shape=exp(k[2]) , scale=k[1]/exp(k[2]) ) , from=0 , to=20 , xlab="lambda (average salamanders)" , ylab="Frequency" )
# same plot as before


################
# ordered logit

# visualize different dists of 4 category data
y <- c(22,26,27,25)
barplot( y , names.arg=1:4 )

y <- c(48,18,16,18)
barplot( y , names.arg=1:4 )

y <- c(16,18,48,18)
barplot( y , names.arg=1:4 )

y <- c(16,18,32,34)
barplot( y , names.arg=1:4 )


# function to compute probability y <= any category value k (1 through max)
# k : collection of values to predict
# z : collection of z intercepts

logit <- function(x) 1/(1+exp(x)) # convert from z to probability
invlogit <- function(x) log(1/x-1) # convert from probability to z

prob.y.k <- function( k , z ) {
    logit( z[k] )
}

y <- c( rep(1,16), rep(2,18), rep(3,32), rep(4,34) ) # fake data: 16 1's, 18 2's, etc.
z <- c( 1 , 0 , -1 , -Inf ) # evenly spread z values
prob.y.k( y , z )

# now show how to fit with mle2
# need a density function, that returns likelihood for each observed value
dorderlogit <- function( x , z , log=FALSE ) {
    p <- logit( z[x] ) # prob y <= k
    nz <- c( Inf , z )
    np <- logit( nz[x] ) # prob y <= k-1
    p <- p - np # subtract to get likelihood of y==k
    if ( log==TRUE ) p <- log(p)
    p
}

# some data to estimate parameters from
y <- c( rep(1,22) , rep(2,26) , rep(3,27) , rep(4,25) )

# perform the fit
m <- mle2( y ~ dorderlogit( z=c(z1,z2,z3,-Inf) ) , start=list( z1=1 , z2=0 , z3=-1 ) )

# plotting the results
plot.new()
frame()
lines( c(0,1) , c(.5,.5) )
lines( c( logit(coef(m)[1]) , logit(coef(m)[1]) ) , c(0,1) ) # z1
lines( c( logit(coef(m)[2]) , logit(coef(m)[2]) ) , c(0,1) ) # z2
lines( c( logit(coef(m)[3]) , logit(coef(m)[3]) ) , c(0,1) ) # z3
lines( c( logit(-Inf) , logit(-Inf) ) , c(0,1) ) # z4 = -Inf

# adding covariates to ordered logit

# covariate values (female/male dummy variable)
f <- c( rep(1,30) , rep(0,30) )
# generate data from a model that makes females have higher average values
y <- rbinom( 60 , prob=1/(1+exp( -f )) , size=4 )
y <- ifelse( y==0 , 1 , y )

# compare female and male means
mean( y[f==1] )
mean( y[f==0] )

# we suppose now that we can add a linear distortion to the z values, so that covariates---like gender---can affect the real cut-offs
z <- c( 1 , 0 , -1 , -Inf )
logit( z ) # male probs
logit( z - 1 ) # female probs --- effect is to increase mean female value

# new density function that allows for linear model distorting sizes of slices
dorderlogit2 <- function( x , z , model , log=FALSE ) {
    p <- logit( z[x] + model ) # prob y <= k
    nz <- c( Inf , z )
    np <- logit( nz[x] + model ) # prob y <= k-1
    p <- p - np # subtract to get likelihood of y==k
    if ( log==TRUE ) p <- log(p)
    p
}

# some data to estimate parameters from
# covariate values (female/male dummy variable)
f <- c( rep(1,30) , rep(0,30) )
# generate data from a model that makes females have higher average values
y <- rbinom( 60 , prob=1/(1+exp( -f )) , size=4 )
y <- ifelse( y==0 , 1 , y )

# perform the fit
m2 <- mle2( y ~ dorderlogit2( z=c(z1,z2,z3,-Inf) , model=b*f ) , start=list( z1=1 , z2=0 , z3=-1 , b=0 ) )

# plotting the results
par( mfcol=c(1,2) )
# males
frame()
lines( c(0,1) , c(.5,.5) )
lines( c( logit(coef(m2)[1]) , logit(coef(m2)[1]) ) , c(0,1) ) # z1
lines( c( logit(coef(m2)[2]) , logit(coef(m2)[2]) ) , c(0,1) ) # z2
lines( c( logit(coef(m2)[3]) , logit(coef(m2)[3]) ) , c(0,1) ) # z3
lines( c( logit(-Inf) , logit(-Inf) ) , c(0,1) ) # z4 = -Inf
# females
frame()
lines( c(0,1) , c(.5,.5) )
lines( c( logit(coef(m2)[1]+coef(m2)[4]) , logit(coef(m2)[1]+coef(m2)[4]) ) , c(0,1) ) # z1
lines( c( logit(coef(m2)[2]+coef(m2)[4]) , logit(coef(m2)[2]+coef(m2)[4]) ) , c(0,1) ) # z2
lines( c( logit(coef(m2)[3]+coef(m2)[4]) , logit(coef(m2)[3]+coef(m2)[4]) ) , c(0,1) ) # z3
lines( c( logit(-Inf) , logit(-Inf) ) , c(0,1) ) # z4 = -Inf

# another way to vizualize
# horizontal axis is our covariate (gender in this case)
# vertical is probability
# plot() call makes the axes and draws the limits for probability of y <= 1
plot( c(0,1) , c( logit( coef(m2)[1] ) , logit( coef(m2)[1]+coef(m2)[4] ) ) , xlab="Gender (female=1)" , ylab="probability" , ylim=c(0,1) , type="l" )
# now draw limits for y <= 2
lines( c(0,1) , c( logit( coef(m2)[2] ) , logit( coef(m2)[2]+coef(m2)[4] ) ) )
# now draw limits for y <= 3
lines( c(0,1) , c( logit( coef(m2)[3] ) , logit( coef(m2)[3]+coef(m2)[4] ) ) )
# now add some text markers for the categories, just for show
text( 0.25 , 0.1 , "1" )
text( 0.35 , 0.3 , "2" )
text( 0.55 , 0.6 , "3" )
text( 0.75 , 0.9 , "4" )

# to highlight one category:
k <- 1
polygon( c(0,1,1,0) , c( logit(coef(m2)[k]) , logit(coef(m2)[k]+coef(m2)[4]) , logit(coef(m2)[k+1]+coef(m2)[4]) , logit(coef(m2)[k+1]) ) , density=10 , col="red" )

# proportional cumulative odds calculation
z2 <- coef(m2)[2]
b <- coef(m2)[4]
# original odds
o1 <- logit(z2)/(1-logit(z2))
# female odds
o2 <- logit(z2+b)/(1-logit(z2+b))
# show same as exp(-b)
o2/o1 # proportional cumulative odds for category 2
exp(-b)

# equivalent via polr() convenience function
# where we used logit( z + b ) above,
# polr uses logit( -z + b ), but otherwise the same
library(MASS)
mpolr <- polr( as.ordered(y) ~ f )