library(languageR)
library(Design)
data(dative)

## Initial exploration
?dative
str(dative)

## give ourselves an indicator variable for the dependent variable
dative$is.PP <- ifelse(dative$RealizationOfRecipient == "PP", 1, 0)

###### Mixed logit model
###### Full model with random effect of verb (item effects)
## Should we edit the model so that it only includes the effects we
## expect to be significant based on previous analyses?
d.glmer1 <- lmer(is.PP ~
                 AccessOfTheme + 
                 AccessOfRec +
                 LengthOfRecipient +
                 AnimacyOfRec +
                 AnimacyOfTheme +
                 PronomOfTheme +
                 DefinOfTheme +
                 LengthOfTheme +
                 SemanticClass +
                 Modality +
                 (1 | Verb), 
                 data = dative, family = "binomial")

## Plot predicted effects of each input
dev.new()
par(mfrow=c(4,3))
plotLMER.fnc(d.glmer1)
par(mfrow=c(4,3))
plotLMER.fnc(d.glmer1, fun=plogis)
par(mfrow=c(1,1))
## How do these compare to the bootstrapped model with cluster
## replacement?  Should we include any slope terms for particular
## verbs?

## Dealing with underdispersion

## Q:  The scale factor in poisson and binomial GLMMs is
## variance/mean.  Variance of a binomial distribution is np(1-p),
## and should therefore scale with the mean, np.  

d.glmer1.q <- lmer(is.PP ~
                   AccessOfTheme + 
                   AccessOfRec +
                   LengthOfRecipient +
                   AnimacyOfRec +
                   AnimacyOfTheme +
                   PronomOfTheme +
                   DefinOfTheme +
                   LengthOfTheme +
                   SemanticClass +
                   Modality +
                   (1 | Verb), 
                   data = dative, family = "quasibinomial")


###### Mixed logit model with crossed random effects
## Full model with random effects for verbs (item effects) and random
## effects for subjects
d.glmer2 <- lmer(is.PP ~
                 AccessOfTheme + 
                 AccessOfRec +
                 LengthOfRecipient +
                 AnimacyOfRec +
                 AnimacyOfTheme +
                 PronomOfTheme +
                 DefinOfTheme +
                 LengthOfTheme +
                 SemanticClass +
                 (1 | Verb) +
                 (1 | Speaker), 
                 data = dative[!is.na(dative$Speaker),],
                 family = "binomial")


## Plot predicted effects of each input
par(mfrow=c(4,3))
plotLMER.fnc(d.glmer2)
par(mfrow=c(4,3))
plotLMER.fnc(d.glmer2, fun=plogis)
par(mfrow=c(1,1))

## Dealing with underdispersion
d.glmer2.q <- lmer(is.PP ~
                   AccessOfTheme + 
                   AccessOfRec +
                   LengthOfRecipient +
                   AnimacyOfRec +
                   AnimacyOfTheme +
                   PronomOfTheme +
                   DefinOfTheme +
                   LengthOfTheme +
                   SemanticClass +
                   (1 | Verb) +
                   (1 | Speaker), 
                   data = dative[!is.na(dative$Speaker),],
                   family = "quasibinomial")



###### Work in progress...
## xyplot(is.PP ~ LengthOfTheme,
##        dative, 
##        type = "p",
##        col = alpha("blue", .5),
##        ylim = c(0,1),
##        panel = function (x, y, ...) {
##          xs = drop(d.glmer1@X)
##          fixed = as.matrix(fixef(d.glmer1))
##          preds = unlist(dimnames(xs)[2])
##          medians = matrix(rep(apply(xs, 2, median), length(x)), ncol=npred, byrow=TRUE)
##          npred = length(fixed)
##          pred.ind = which(preds == "LengthOfTheme")
##          if (pred.ind == npred) {
##            panel.curve(plogis(cbind(medians[,1:(pred.ind-1)], x) %*% fixed))
##          }
##          else {
##            print(length(x))
##            print(dim(cbind(medians[,1:(pred.ind-1)], x, medians[,(pred.ind+1):npred])))
##            print(dim(fixed))
##            print(dim(cbind(medians[,1:(pred.ind-1)], x, medians[,(pred.ind+1):npred]) %*% fixed))
##            panel.curve(plogis(cbind(medians[,1:(pred.ind-1)], x, medians[,(pred.ind+1):npred]) %*% fixed))
##          }
##          panel.stripplot(x, y, jitter=TRUE, ...)
##        })

###### Simulate the model
## copied from slot(getMethods("simulate"), "methods")$mer
## my.simulate <- function (object, nsim = 1, seed = NULL, ...) 
## {
##   if (!exists(".Random.seed", envir = .GlobalEnv)) 
##     runif(1)
##   if (is.null(seed)) 
##     RNGstate <- .Random.seed
##   else {
##     R.seed <- .Random.seed
##     set.seed(seed)
##     RNGstate <- structure(seed, kind = as.list(RNGkind()))
##     on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
##   }
##   stopifnot((nsim <- as.integer(nsim[1])) > 0, is(object, "lmer"))
##   lpred <- .Call(lme4:::mer_simulate, object, nsim)
##   sc <- abs(object@devComp[8])
##   lpred <- lpred + drop(object@X %*% fixef(object))
##   n <- prod(dim(lpred))
##   fam <- attr(object, "family")
##   if (fam$family == "binomial") {
##     summary(lpred)
##     response <- as.data.frame(matrix(byrow = FALSE, ncol = nsim, 
##                                      rbinom(n, 1, fam$linkinv(lpred))))
##     attr(response, "seed") <- RNGstate
##     return(response)
##   }
## }