Modelos de Regressão para Dados de Contagem com o R

github.com/leg-ufpr/MRDCr

Análise exploratória

library(MRDCr)
# help(soja)
ls("package:MRDCr")
##  [1] "apc"                  "calc_mean_cmp"       
##  [3] "calc_mean_gcnt"       "calc_var_cmp"        
##  [5] "cambras"              "capdesfo"            
##  [7] "capmosca"             "cholV_eta"           
##  [9] "cmp"                  "conftemp"            
## [11] "confterm"             "convergencez"        
## [13] "dcmp"                 "dgcnt"               
## [15] "dpgnz"                "gcnt"                
## [17] "led"                  "llcmp"               
## [19] "llgcnt"               "llpgnz"              
## [21] "nematoide"            "ninfas"              
## [23] "panel.beeswarm"       "panel.cbH"           
## [25] "panel.groups.segplot" "peixe"               
## [27] "pgnz"                 "postura"             
## [29] "predict.mle2"         "prepanel.cbH"        
## [31] "seguros"              "soja"
library(lattice)

str(soja)
## 'data.frame':    75 obs. of  5 variables:
##  $ K   : int  0 30 60 120 180 0 30 60 120 180 ...
##  $ umid: Factor w/ 3 levels "37,5","50","62,5": 1 1 1 1 1 2 2 2 2 2 ...
##  $ bloc: Factor w/ 5 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ ngra: int  136 159 156 171 190 140 193 200 208 237 ...
##  $ nvag: int  56 62 66 68 82 63 86 94 86 97 ...
xtabs(~umid + K, data = soja[-75, ])
##       K
## umid   0 30 60 120 180
##   37,5 5  5  5   5   5
##   50   5  5  5   5   5
##   62,5 5  5  5   5   4
xyplot(nvag + ngra ~ K, groups = umid, outer = TRUE,
       data = soja[-74, ],
       type = c("p", "a"), scales = "free",
       ylab = NULL,
       xlab = expression("Nível de potássio aplicado"~(mg~dm^{-3})),
       auto.key = list(title = "Umidade do solo (%)",
                       cex.title = 1,
                       columns = 3),
       strip = strip.custom(
           factor.levels = c("Número de vagens",
                             "Número de grãos")))

soja <- soja[-74, ]

GLM Poisson para número de vagens viáveis por vaso

# Considerar K categórico.
soja <- transform(soja, K = factor(K))

m0 <- glm(nvag ~ bloc + umid * K, data = soja, family = poisson)

par(mfrow = c(2, 2))
plot(m0); layout(1)

#-----------------------------------------------------------------------

deviance(m0)
## [1] 65.27834
df.residual(m0)
## [1] 55
summary(m0)
## 
## Call:
## glm(formula = nvag ~ bloc + umid * K, family = poisson, data = soja)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9681  -0.6393  -0.0195   0.4854   2.3306  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.95369    0.06892  57.363  < 2e-16 ***
## blocII        -0.02928    0.04091  -0.716  0.47414    
## blocIII       -0.07265    0.04136  -1.756  0.07902 .  
## blocIV        -0.12544    0.04194  -2.991  0.00278 ** 
## blocV         -0.10795    0.04297  -2.512  0.01199 *  
## umid50         0.13404    0.08765   1.529  0.12619    
## umid62,5       0.21656    0.08602   2.518  0.01181 *  
## K30            0.27427    0.08493   3.229  0.00124 ** 
## K60            0.30797    0.08432   3.652  0.00026 ***
## K120           0.32883    0.08395   3.917 8.97e-05 ***
## K180           0.25540    0.08528   2.995  0.00275 ** 
## umid50:K30     0.06322    0.11557   0.547  0.58433    
## umid62,5:K30  -0.10747    0.11536  -0.932  0.35152    
## umid50:K60     0.16561    0.11370   1.457  0.14521    
## umid62,5:K60   0.10735    0.11220   0.957  0.33869    
## umid50:K120    0.14920    0.11338   1.316  0.18818    
## umid62,5:K120  0.11839    0.11404   1.038  0.29922    
## umid50:K180    0.30370    0.11361   2.673  0.00751 ** 
## umid62,5:K180  0.19838    0.11256   1.762  0.07800 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 322.684  on 73  degrees of freedom
## Residual deviance:  65.278  on 55  degrees of freedom
## AIC: 557.23
## 
## Number of Fisher Scoring iterations: 4
anova(m0, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: nvag
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                      73     322.68              
## bloc    4   14.284        69     308.40  0.006442 ** 
## umid    2   92.911        67     215.49 < 2.2e-16 ***
## K       4  136.060        63      79.43 < 2.2e-16 ***
## umid:K  8   14.150        55      65.28  0.077926 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Predição com intervalos de confiança.

library(doBy)
library(multcomp)

X <- LSmatrix(m0, effect = c("umid", "K"))

pred <- attr(X, "grid")
pred <- transform(pred,
                  K = as.integer(as.character(K)),
                  umid = factor(umid))

# Quantil normal para fazer um IC de 95%.
qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)

# Preditos pela Poisson.
# aux <- predict(m0, newdata = pred$pois, se.fit = TRUE)
# aux <- exp(aux$fit + outer(aux$se.fit, qn, FUN = "*"))
# pred$pois <- cbind(pred$pois, aux)
aux <- confint(glht(m0, linfct = X),
               calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred <- cbind(pred, exp(aux))
str(pred)
## 'data.frame':    15 obs. of  5 variables:
##  $ umid: Factor w/ 3 levels "37,5","50","62,5": 1 2 3 1 2 3 1 2 3 1 ...
##  $ K   : int  0 0 0 30 30 30 60 60 60 120 ...
##  $ fit : num  48.7 55.7 60.5 64.1 78.1 ...
##  $ lwr : num  43 49.6 54.1 57.5 70.7 ...
##  $ upr : num  55.3 62.7 67.7 71.5 86.3 ...
xyplot(fit ~ K | umid, data = pred,
       layout = c(NA, 1), as.table = TRUE,
       ylab = "Número de vagens por vaso",
       xlab = expression("Nível de potássio aplicado"~(mg~dm^{-3})),
       xlim = extendrange(range(pred$K), f = 0.1),
       ly = pred$lwr, uy = pred$upr, cty = "bars", length = 0,
       prepanel = prepanel.cbH,
       panel = panel.cbH)

#-----------------------------------------------------------------------
# Comparações múltiplas.

L <- by(X, INDICES = pred$umid, FUN = as.matrix)
names(L) <- levels(soja$umid)
L <- lapply(L, "rownames<-", levels(soja$K))
K <- lapply(L, apc)
apc(L[[1]])
##         (Intercept) blocII blocIII blocIV blocV umid50 umid62,5 K30
## 0-30              0      0       0      0     0      0        0  -1
## 0-60              0      0       0      0     0      0        0   0
## 0-120             0      0       0      0     0      0        0   0
## 0-180             0      0       0      0     0      0        0   0
## 30-60             0      0       0      0     0      0        0   1
## 30-120            0      0       0      0     0      0        0   1
## 30-180            0      0       0      0     0      0        0   1
## 60-120            0      0       0      0     0      0        0   0
## 60-180            0      0       0      0     0      0        0   0
## 120-180           0      0       0      0     0      0        0   0
##         K60 K120 K180 umid50:K30 umid62,5:K30 umid50:K60
## 0-30      0    0    0          0            0          0
## 0-60     -1    0    0          0            0          0
## 0-120     0   -1    0          0            0          0
## 0-180     0    0   -1          0            0          0
## 30-60    -1    0    0          0            0          0
## 30-120    0   -1    0          0            0          0
## 30-180    0    0   -1          0            0          0
## 60-120    1   -1    0          0            0          0
## 60-180    1    0   -1          0            0          0
## 120-180   0    1   -1          0            0          0
##         umid62,5:K60 umid50:K120 umid62,5:K120 umid50:K180
## 0-30               0           0             0           0
## 0-60               0           0             0           0
## 0-120              0           0             0           0
## 0-180              0           0             0           0
## 30-60              0           0             0           0
## 30-120             0           0             0           0
## 30-180             0           0             0           0
## 60-120             0           0             0           0
## 60-180             0           0             0           0
## 120-180            0           0             0           0
##         umid62,5:K180
## 0-30                0
## 0-60                0
## 0-120               0
## 0-180               0
## 30-60               0
## 30-120              0
## 30-180              0
## 60-120              0
## 60-180              0
## 120-180             0
lapply(K,
       FUN = function(k) {
           summary(glht(model = m0, linfct = k),
                   test = adjusted(type = "fdr"))
       })
## $`37,5`
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = nvag ~ bloc + umid * K, family = poisson, data = soja)
## 
## Linear Hypotheses:
##              Estimate Std. Error z value Pr(>|z|)    
## 0-30 == 0    -0.27427    0.08493  -3.229 0.004137 ** 
## 0-60 == 0    -0.30797    0.08432  -3.652 0.001300 ** 
## 0-120 == 0   -0.32883    0.08395  -3.917 0.000897 ***
## 0-180 == 0   -0.25540    0.08528  -2.995 0.006865 ** 
## 30-60 == 0   -0.03369    0.07828  -0.430 0.811949    
## 30-120 == 0  -0.05456    0.07788  -0.701 0.719952    
## 30-180 == 0   0.01887    0.07931   0.238 0.811949    
## 60-120 == 0  -0.02087    0.07721  -0.270 0.811949    
## 60-180 == 0   0.05256    0.07866   0.668 0.719952    
## 120-180 == 0  0.07343    0.07826   0.938 0.696218    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
## 
## 
## $`50`
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = nvag ~ bloc + umid * K, family = poisson, data = soja)
## 
## Linear Hypotheses:
##               Estimate Std. Error z value Pr(>|z|)    
## 0-30 == 0    -0.337496   0.078369  -4.306 4.15e-05 ***
## 0-60 == 0    -0.473581   0.076265  -6.210 1.77e-09 ***
## 0-120 == 0   -0.478036   0.076200  -6.273 1.77e-09 ***
## 0-180 == 0   -0.559104   0.075056  -7.449 9.39e-13 ***
## 30-60 == 0   -0.136086   0.069208  -1.966  0.07037 .  
## 30-120 == 0  -0.140540   0.069136  -2.033  0.07012 .  
## 30-180 == 0  -0.221608   0.067873  -3.265  0.00219 ** 
## 60-120 == 0  -0.004454   0.066741  -0.067  0.94679    
## 60-180 == 0  -0.085522   0.065432  -1.307  0.23870    
## 120-180 == 0 -0.081068   0.065356  -1.240  0.23870    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
## 
## 
## $`62,5`
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = nvag ~ bloc + umid * K, family = poisson, data = soja)
## 
## Linear Hypotheses:
##               Estimate Std. Error z value Pr(>|z|)    
## 0-30 == 0    -0.166800   0.078062  -2.137 0.046595 *  
## 0-60 == 0    -0.415317   0.074020  -5.611 6.71e-08 ***
## 0-120 == 0   -0.447218   0.077181  -5.794 3.43e-08 ***
## 0-180 == 0   -0.453784   0.073463  -6.177 6.53e-09 ***
## 30-60 == 0   -0.248517   0.070512  -3.524 0.000707 ***
## 30-120 == 0  -0.280417   0.073823  -3.799 0.000291 ***
## 30-180 == 0  -0.286984   0.069927  -4.104 0.000101 ***
## 60-120 == 0  -0.031900   0.069536  -0.459 0.718229    
## 60-180 == 0  -0.038466   0.065384  -0.588 0.695404    
## 120-180 == 0 -0.006566   0.068942  -0.095 0.924124    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)

GLM Poisson para número de grãos por vaso

m1 <- glm(ngra ~ bloc + umid * K, data = soja, family = poisson)

par(mfrow = c(2, 2))
plot(m1); layout(1)

#-----------------------------------------------------------------------

deviance(m1)
## [1] 124.7194
df.residual(m1)
## [1] 55
summary(m1)
## 
## Call:
## glm(formula = ngra ~ bloc + umid * K, family = poisson, data = soja)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6297  -1.0707  -0.0873   0.7428   3.2810  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    4.80196    0.04477 107.254  < 2e-16 ***
## blocII        -0.01939    0.02655  -0.730 0.465260    
## blocIII       -0.03663    0.02667  -1.373 0.169673    
## blocIV        -0.10559    0.02715  -3.889 0.000101 ***
## blocV         -0.09313    0.02788  -3.340 0.000837 ***
## umid50         0.13245    0.05692   2.327 0.019968 *  
## umid62,5       0.18548    0.05623   3.299 0.000972 ***
## K30            0.29799    0.05486   5.432 5.56e-08 ***
## K60            0.34434    0.05432   6.339 2.32e-10 ***
## K120           0.36493    0.05409   6.746 1.52e-11 ***
## K180           0.29542    0.05489   5.383 7.35e-08 ***
## umid50:K30     0.04344    0.07482   0.581 0.561495    
## umid62,5:K30  -0.13669    0.07527  -1.816 0.069349 .  
## umid50:K60     0.11559    0.07361   1.570 0.116354    
## umid62,5:K60   0.09174    0.07289   1.259 0.208190    
## umid50:K120    0.11860    0.07329   1.618 0.105637    
## umid62,5:K120  0.16266    0.07367   2.208 0.027242 *  
## umid50:K180    0.28832    0.07327   3.935 8.33e-05 ***
## umid62,5:K180  0.21569    0.07285   2.961 0.003071 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 761.43  on 73  degrees of freedom
## Residual deviance: 124.72  on 55  degrees of freedom
## AIC: 681.34
## 
## Number of Fisher Scoring iterations: 4
anova(m1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: ngra
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                      73     761.43              
## bloc    4    28.34        69     733.09 1.063e-05 ***
## umid    2   185.06        67     548.03 < 2.2e-16 ***
## K       4   380.32        63     167.71 < 2.2e-16 ***
## umid:K  8    42.99        55     124.72 8.830e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Predição com intervalos de confiança.

X <- LSmatrix(m1, effect = c("umid", "K"))

pred <- attr(X, "grid")
pred <- transform(pred,
                  K = as.integer(as.character(K)),
                  umid = factor(umid))

# Quantil normal para fazer um IC de 95%.
qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)

aux <- confint(glht(m1, linfct = X),
               calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred <- cbind(pred, exp(aux))
str(pred)
## 'data.frame':    15 obs. of  5 variables:
##  $ umid: Factor w/ 3 levels "37,5","50","62,5": 1 2 3 1 2 3 1 2 3 1 ...
##  $ K   : int  0 0 0 30 30 30 60 60 60 120 ...
##  $ fit : num  116 132 139 156 186 ...
##  $ lwr : num  107 122 129 145 174 ...
##  $ upr : num  126 143 150 167 198 ...
xyplot(fit ~ K | umid, data = pred,
       layout = c(NA, 1), as.table = TRUE,
       ylab = "Número de grãos por vaso",
       xlab = expression("Nível de potássio aplicado"~(mg~dm^{-3})),
       xlim = extendrange(range(pred$K), f = 0.1),
       ly = pred$lwr, uy = pred$upr, cty = "bars", length = 0,
       prepanel = prepanel.cbH,
       panel = panel.cbH)

#-----------------------------------------------------------------------
# Comparações múltiplas.

L <- by(X, INDICES = pred$umid, FUN = as.matrix)
names(L) <- levels(soja$umid)
L <- lapply(L, "rownames<-", levels(soja$K))
K <- lapply(L, apc)
apc(L[[1]])
##         (Intercept) blocII blocIII blocIV blocV umid50 umid62,5 K30
## 0-30              0      0       0      0     0      0        0  -1
## 0-60              0      0       0      0     0      0        0   0
## 0-120             0      0       0      0     0      0        0   0
## 0-180             0      0       0      0     0      0        0   0
## 30-60             0      0       0      0     0      0        0   1
## 30-120            0      0       0      0     0      0        0   1
## 30-180            0      0       0      0     0      0        0   1
## 60-120            0      0       0      0     0      0        0   0
## 60-180            0      0       0      0     0      0        0   0
## 120-180           0      0       0      0     0      0        0   0
##         K60 K120 K180 umid50:K30 umid62,5:K30 umid50:K60
## 0-30      0    0    0          0            0          0
## 0-60     -1    0    0          0            0          0
## 0-120     0   -1    0          0            0          0
## 0-180     0    0   -1          0            0          0
## 30-60    -1    0    0          0            0          0
## 30-120    0   -1    0          0            0          0
## 30-180    0    0   -1          0            0          0
## 60-120    1   -1    0          0            0          0
## 60-180    1    0   -1          0            0          0
## 120-180   0    1   -1          0            0          0
##         umid62,5:K60 umid50:K120 umid62,5:K120 umid50:K180
## 0-30               0           0             0           0
## 0-60               0           0             0           0
## 0-120              0           0             0           0
## 0-180              0           0             0           0
## 30-60              0           0             0           0
## 30-120             0           0             0           0
## 30-180             0           0             0           0
## 60-120             0           0             0           0
## 60-180             0           0             0           0
## 120-180            0           0             0           0
##         umid62,5:K180
## 0-30                0
## 0-60                0
## 0-120               0
## 0-180               0
## 30-60               0
## 30-120              0
## 30-180              0
## 60-120              0
## 60-180              0
## 120-180             0
lapply(K,
       FUN = function(k) {
           summary(glht(model = m1, linfct = k),
                   test = adjusted(type = "fdr"))
       })
## $`37,5`
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = ngra ~ bloc + umid * K, family = poisson, data = soja)
## 
## Linear Hypotheses:
##               Estimate Std. Error z value Pr(>|z|)    
## 0-30 == 0    -0.297991   0.054856  -5.432 1.84e-07 ***
## 0-60 == 0    -0.344337   0.054324  -6.339 1.16e-09 ***
## 0-120 == 0   -0.364931   0.054094  -6.746 1.52e-10 ***
## 0-180 == 0   -0.295424   0.054886  -5.383 1.84e-07 ***
## 30-60 == 0   -0.046345   0.050060  -0.926    0.443    
## 30-120 == 0  -0.066939   0.049811  -1.344    0.298    
## 30-180 == 0   0.002567   0.050670   0.051    0.960    
## 60-120 == 0  -0.020594   0.049224  -0.418    0.751    
## 60-180 == 0   0.048913   0.050093   0.976    0.443    
## 120-180 == 0  0.069507   0.049844   1.394    0.298    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
## 
## 
## $`50`
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = ngra ~ bloc + umid * K, family = poisson, data = soja)
## 
## Linear Hypotheses:
##              Estimate Std. Error z value Pr(>|z|)    
## 0-30 == 0    -0.34143    0.05087  -6.711 4.82e-11 ***
## 0-60 == 0    -0.45993    0.04968  -9.258  < 2e-16 ***
## 0-120 == 0   -0.48353    0.04945  -9.777  < 2e-16 ***
## 0-180 == 0   -0.58374    0.04855 -12.024  < 2e-16 ***
## 30-60 == 0   -0.11850    0.04506  -2.630  0.01068 *  
## 30-120 == 0  -0.14210    0.04481  -3.171  0.00253 ** 
## 30-180 == 0  -0.24231    0.04381  -5.531 6.36e-08 ***
## 60-120 == 0  -0.02360    0.04345  -0.543  0.58707    
## 60-180 == 0  -0.12381    0.04241  -2.919  0.00501 ** 
## 120-180 == 0 -0.10022    0.04215  -2.378  0.01936 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
## 
## 
## $`62,5`
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glm(formula = ngra ~ bloc + umid * K, family = poisson, data = soja)
## 
## Linear Hypotheses:
##              Estimate Std. Error z value Pr(>|z|)    
## 0-30 == 0    -0.16130    0.05153  -3.130   0.0025 ** 
## 0-60 == 0    -0.43608    0.04860  -8.972  < 2e-16 ***
## 0-120 == 0   -0.52759    0.05001 -10.550  < 2e-16 ***
## 0-180 == 0   -0.51111    0.04791 -10.668  < 2e-16 ***
## 30-60 == 0   -0.27478    0.04635  -5.928 5.11e-09 ***
## 30-120 == 0  -0.36629    0.04782  -7.659 3.73e-14 ***
## 30-180 == 0  -0.34981    0.04562  -7.667 3.73e-14 ***
## 60-120 == 0  -0.09152    0.04465  -2.050   0.0505 .  
## 60-180 == 0  -0.07504    0.04229  -1.774   0.0844 .  
## 120-180 == 0  0.01648    0.04389   0.375   0.7073    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)

Informações da sessão

## Atualizado em 21 de julho de 2016.
## 
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04 LTS
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=pt_BR.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets 
## [7] methods   base     
## 
## other attached packages:
##  [1] multcomp_1.4-6      TH.data_1.0-7       MASS_7.3-45        
##  [4] survival_2.39-5     mvtnorm_1.0-5       doBy_4.5-15        
##  [7] rmarkdown_1.0       latticeExtra_0.6-28 RColorBrewer_1.1-2 
## [10] lattice_0.20-33     knitr_1.13          MRDCr_0.0-2        
## [13] gsubfn_0.6-6        proto_0.3-10        bbmle_1.0.18       
## [16] devtools_1.12.0    
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.5       formatR_1.4       compiler_3.3.1   
##  [4] tools_3.3.1       testthat_1.0.2    digest_0.6.9     
##  [7] evaluate_0.9      memoise_1.0.0     Matrix_1.2-6     
## [10] yaml_2.1.13       withr_1.0.2       stringr_1.0.0    
## [13] roxygen2_5.0.1    grid_3.3.1        R6_2.1.2         
## [16] tcltk_3.3.1       magrittr_1.5      htmltools_0.3.5  
## [19] codetools_0.2-14  splines_3.3.1     numDeriv_2014.2-1
## [22] sandwich_2.3-4    stringi_1.1.1     crayon_1.3.2     
## [25] zoo_1.7-13