1 library(languageR)
   2 library(Design)
   3 data(dative)
   4 
   5 ## Initial exploration
   6 ?dative
   7 str(dative)
   8 
   9 ## give ourselves an indicator variable for the dependent variable
  10 dative$is.PP <- ifelse(dative$RealizationOfRecipient == "PP", 1, 0)
  11 
  12 ###### Mixed logit model
  13 ###### Full model with random effect of verb (item effects)
  14 ## Should we edit the model so that it only includes the effects we
  15 ## expect to be significant based on previous analyses?
  16 d.glmer1 <- lmer(is.PP ~
  17                  AccessOfTheme + 
  18                  AccessOfRec +
  19                  LengthOfRecipient +
  20                  AnimacyOfRec +
  21                  AnimacyOfTheme +
  22                  PronomOfTheme +
  23                  DefinOfTheme +
  24                  LengthOfTheme +
  25                  SemanticClass +
  26                  Modality +
  27                  (1 | Verb), 
  28                  data = dative, family = "binomial")
  29 
  30 ## Plot predicted effects of each input
  31 dev.new()
  32 par(mfrow=c(4,3))
  33 plotLMER.fnc(d.glmer1)
  34 par(mfrow=c(4,3))
  35 plotLMER.fnc(d.glmer1, fun=plogis)
  36 par(mfrow=c(1,1))
  37 ## How do these compare to the bootstrapped model with cluster
  38 ## replacement?  Should we include any slope terms for particular
  39 ## verbs?
  40 
  41 ## Dealing with underdispersion
  42 
  43 ## Q:  The scale factor in poisson and binomial GLMMs is
  44 ## variance/mean.  Variance of a binomial distribution is np(1-p),
  45 ## and should therefore scale with the mean, np.  
  46 
  47 d.glmer1.q <- lmer(is.PP ~
  48                    AccessOfTheme + 
  49                    AccessOfRec +
  50                    LengthOfRecipient +
  51                    AnimacyOfRec +
  52                    AnimacyOfTheme +
  53                    PronomOfTheme +
  54                    DefinOfTheme +
  55                    LengthOfTheme +
  56                    SemanticClass +
  57                    Modality +
  58                    (1 | Verb), 
  59                    data = dative, family = "quasibinomial")
  60 
  61 
  62 ###### Mixed logit model with crossed random effects
  63 ## Full model with random effects for verbs (item effects) and random
  64 ## effects for subjects
  65 d.glmer2 <- lmer(is.PP ~
  66                  AccessOfTheme + 
  67                  AccessOfRec +
  68                  LengthOfRecipient +
  69                  AnimacyOfRec +
  70                  AnimacyOfTheme +
  71                  PronomOfTheme +
  72                  DefinOfTheme +
  73                  LengthOfTheme +
  74                  SemanticClass +
  75                  (1 | Verb) +
  76                  (1 | Speaker), 
  77                  data = dative[!is.na(dative$Speaker),],
  78                  family = "binomial")
  79 
  80 
  81 ## Plot predicted effects of each input
  82 par(mfrow=c(4,3))
  83 plotLMER.fnc(d.glmer2)
  84 par(mfrow=c(4,3))
  85 plotLMER.fnc(d.glmer2, fun=plogis)
  86 par(mfrow=c(1,1))
  87 
  88 ## Dealing with underdispersion
  89 d.glmer2.q <- lmer(is.PP ~
  90                    AccessOfTheme + 
  91                    AccessOfRec +
  92                    LengthOfRecipient +
  93                    AnimacyOfRec +
  94                    AnimacyOfTheme +
  95                    PronomOfTheme +
  96                    DefinOfTheme +
  97                    LengthOfTheme +
  98                    SemanticClass +
  99                    (1 | Verb) +
 100                    (1 | Speaker), 
 101                    data = dative[!is.na(dative$Speaker),],
 102                    family = "quasibinomial")
 103 
 104 
 105 
 106 ###### Work in progress...
 107 ## xyplot(is.PP ~ LengthOfTheme,
 108 ##        dative, 
 109 ##        type = "p",
 110 ##        col = alpha("blue", .5),
 111 ##        ylim = c(0,1),
 112 ##        panel = function (x, y, ...) {
 113 ##          xs = drop(d.glmer1@X)
 114 ##          fixed = as.matrix(fixef(d.glmer1))
 115 ##          preds = unlist(dimnames(xs)[2])
 116 ##          medians = matrix(rep(apply(xs, 2, median), length(x)), ncol=npred, byrow=TRUE)
 117 ##          npred = length(fixed)
 118 ##          pred.ind = which(preds == "LengthOfTheme")
 119 ##          if (pred.ind == npred) {
 120 ##            panel.curve(plogis(cbind(medians[,1:(pred.ind-1)], x) %*% fixed))
 121 ##          }
 122 ##          else {
 123 ##            print(length(x))
 124 ##            print(dim(cbind(medians[,1:(pred.ind-1)], x, medians[,(pred.ind+1):npred])))
 125 ##            print(dim(fixed))
 126 ##            print(dim(cbind(medians[,1:(pred.ind-1)], x, medians[,(pred.ind+1):npred]) %*% fixed))
 127 ##            panel.curve(plogis(cbind(medians[,1:(pred.ind-1)], x, medians[,(pred.ind+1):npred]) %*% fixed))
 128 ##          }
 129 ##          panel.stripplot(x, y, jitter=TRUE, ...)
 130 ##        })
 131 
 132 ###### Simulate the model
 133 ## copied from slot(getMethods("simulate"), "methods")$mer
 134 ## my.simulate <- function (object, nsim = 1, seed = NULL, ...) 
 135 ## {
 136 ##   if (!exists(".Random.seed", envir = .GlobalEnv)) 
 137 ##     runif(1)
 138 ##   if (is.null(seed)) 
 139 ##     RNGstate <- .Random.seed
 140 ##   else {
 141 ##     R.seed <- .Random.seed
 142 ##     set.seed(seed)
 143 ##     RNGstate <- structure(seed, kind = as.list(RNGkind()))
 144 ##     on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
 145 ##   }
 146 ##   stopifnot((nsim <- as.integer(nsim[1])) > 0, is(object, "lmer"))
 147 ##   lpred <- .Call(lme4:::mer_simulate, object, nsim)
 148 ##   sc <- abs(object@devComp[8])
 149 ##   lpred <- lpred + drop(object@X %*% fixef(object))
 150 ##   n <- prod(dim(lpred))
 151 ##   fam <- attr(object, "family")
 152 ##   if (fam$family == "binomial") {
 153 ##     summary(lpred)
 154 ##     response <- as.data.frame(matrix(byrow = FALSE, ncol = nsim, 
 155 ##                                      rbinom(n, 1, fam$linkinv(lpred))))
 156 ##     attr(response, "seed") <- RNGstate
 157 ##     return(response)
 158 ##   }
 159 ## }