Movatterモバイル変換

A Logit-Normal GLMM Dataset

Description

This data set contains simulated data from the paper of Booth and Hobert (referenced below) as well as another vector.

Usage

data(Booth2)

Format

A data frame with 3 columns:

y

Response vector.

x1

Fixed effect model matrix. The matrix has just one column vector.

z1

A categorical vector to be used for part of the random effect model matrix.

z2

A categorical vector to be used for part of the random effect model matrix.

Details

The original data set was generated by Booth and Hobert using a single variance component, a single fixed effect, no intercept, and a logit link. This data set has the z2 vector added purely to illustrate an example with multiple variance components.

References

Booth, J. G. and Hobert, J. P. (1999) Maximizing generalized linear mixed model likelihoods with an automated Monte Carlo EM algorithm.Journal of the Royal Statistical Society, Series B,61, 265–285.doi:10.1111/1467-9868.00176.

Examples

data(Booth2)

A Logit-Normal GLMM Dataset from Booth and Hobert

Description

This data set contains simulated data from the paper of Booth and Hobert referenced below.

Usage

data(BoothHobert)

Format

A data frame with 3 columns:

y

Response vector.

x1

Fixed effect model matrix. The matrix has just one column vector.

z1

Random effect model matrix. The matrix has just one column vector.

Details

This data set was generated by Booth and Hobert using a single variance component, a single fixed effect, no intercept, and a logit link.

References

Booth, J. G. and Hobert, J. P. (1999) Maximizing generalized linear mixed model likelihoods with an automated Monte Carlo EM algorithm.Journal of the Royal Statistical Society, Series B,61, 265–285.doi:10.1111/1467-9868.00176.

Examples

data(BoothHobert)

Presence of Bacteria after Drug Treatments

Description

Tests of the presence of the bacteria H. influenzae in children with otitis media in the Northern Territory of Australia.

Usage

data(bacteria)

Format

A data frame with the following columns:

y

Presence or absence: a factor with levelsn andy.

ap

active/placebo: a factor with levelsa andp.

hilo

hi/low compliance: a factor with levelshi andlo.

week

Numeric: week of test.

ID

Subject ID: a factor.

trt

A factor with levelsplacebo,drug,drug+, a re-coding ofap andhilo.

y2

y reformatted as 0/1 rather than n/y.

Details

Dr. A. Leach tested the effects of a drug on 50 children with a history of otitis media in the Northern Territory of Australia. The children were randomized to the drug or the a placebo, and also to receive active encouragement to comply with taking the drug.

The presence of H. influenzae was checked at weeks 0, 2, 4, 6 and 11: 30 of the checks were missing and are not included in this data frame.

References

Venables, W. N. and Ripley, B. D. (2002)Modern Applied Statistics with S, Fourth edition. Springer.

Examples

data(bacteria)

Functions for the Bernoulli family.

Description

Given a scalareta, this calculates the cumulant and two derivatives for the Bernoulli family. Also checks that the data are entered correctly.

Usage

bernoulli.glmm()

Value

Note

This function is to be used by theglmm command.

Author(s)

Christina Knudson

See Also

glmm

Examples

eta<--3:3bernoulli.glmm()$family.glmmbernoulli.glmm()$cum(eta)bernoulli.glmm()$cp(1)bernoulli.glmm()$cpp(2)

Functions for the Binomial family.

Description

Given a scalareta and the number of trials, this calculates the cumulant and two derivatives for the Bernoulli family. Also checks that the data are entered correctly.

Usage

binomial.glmm()

Value

Note

This function is to be used by theglmm command.

Author(s)

Christina Knudson

See Also

glmm

Examples

eta<--3:3ntrials <- 1binomial.glmm()$family.glmmbinomial.glmm()$cum(eta, ntrials)binomial.glmm()$cp(1, ntrials)binomial.glmm()$cpp(2, ntrials)

Contagious Bovine Pleuropneumonia

Description

This data set is a reformatted version ofcbpp from thelme4 package. Contagious bovine pleuropneumonia (CBPP) is a major disease of cattle in Africa, caused by a mycoplasma. This dataset describes the serological incidence of CBPP in zebu cattle during a follow-up survey implemented in 15 commercial herds located in the Boji district of Ethiopia. The goal of the survey was to study the within-herd spread of CBPP in newly infected herds. Blood samples were quarterly collected from all animals of these herds to determine their CBPP status. These data were used to compute the serological incidence of CBPP (new cases occurring during a given time period). Some data are missing (lost to follow-up).

Usage

data(cbpp2)

Format

A data frame with 3 columns:

Y

Response vector. 1 if CBPP is observed, 0 otherwise.

period

A factor with levels 1 to 4.

herd

A factor identifying the herd (1 through 15).

Details

Serological status was determined using a competitive enzyme-linked immuno-sorbent assay (cELISA).

References

Lesnoff, M., Laval, G., Bonnet, P., Abdicho, S., Workalemahu, A., Kifle, D., Peyraud, A., Lancelot, R., Thiaucourt, F. (2004) Within-herd spread of contagious bovine pleuropneumonia in Ethiopian highlands.Preventive Veterinary Medicine,64, 27–40.doi:10.1016/j.prevetmed.2004.03.005.

Examples

data(cbpp2)

Extract Model Coefficients

Description

A function that extracts the fixed effect coefficients returned fromglmm.

Usage

## S3 method for class 'glmm'coef(object,...)

Arguments

Value

Author(s)

Christina Knudson

See Also

glmm for model fitting.

Examples

library(glmm)set.seed(1234)data(salamander)m<-1000sal<-glmm(Mate~0+Cross,random=list(~0+Female,~0+Male),varcomps.names=c("F","M"),data=salamander,family.glmm=bernoulli.glmm,m=m,debug=TRUE,doPQL=FALSE)coef(sal)

Calculates Asymptotic Confidence Intervals

Description

A function that calculates asymptotic confidence intervals for one or more parameters in a model fitted by byglmm. Confidence intervals can be calculated for fixed effect parameters and variance components using models.

Usage

## S3 method for class 'glmm'confint(object, parm, level, ...)

Arguments

Value

A matrix (or vector) with columns giving lower and upper confidence limits for each parameter. These will be labeled as (1-level)/2 and 1-(1-level)/2 in percent. By default, 2.5

Author(s)

Christina Knudson

See Also

glmm for model fitting.

Examples

library(glmm)data(BoothHobert)set.seed(123)mod.mcml1<-glmm(y~0+x1,list(y~0+z1),varcomps.names=c("z1"), data=BoothHobert,family.glmm=bernoulli.glmm,m=1000,doPQL=TRUE)confint(mod.mcml1)

Fitting Generalized Linear Mixed Models using MCML

Description

This function fits generalized linear mixed models (GLMMs) by approximating the likelihood with ordinary Monte Carlo, then maximizing the approximated likelihood.

Usage

glmm(fixed, random, varcomps.names, data, family.glmm, m, varcomps.equal, weights=NULL, doPQL = TRUE,debug=FALSE, p1=1/3,p2=1/3, p3=1/3,rmax=1000,iterlim=1000, par.init, zeta=5, cluster=NULL)

Arguments

Details

Let\beta be a vector of fixed effects and letu be a vector of random effects. LetX andZ be design matrices for the fixed and random effects, respectively. The random effects are assumed to be normally distributed with mean 0 and variance matrixD, whereD is diagonal with entries from the unknown vector\nu. Lettingg be the link function,g(\mu)=X \beta+ZU. If the response type is Bernoulli or Binomial, then the logit function is the link; if the response type is Poisson, then the natural logarithm is the link function.

Models for glmm are specified symbolically. A typical fixed effects model has the formresponse ~ terms whereresponse is the (numeric) response vector andterms is a series of terms which specifies a linear predictor for response. A terms specification of the formfirst + second indicates all the terms infirst together with all the terms insecond with duplicates removed.

A specification of the formfirst:second indicates the set of terms obtained by taking the interactions of all terms infirst with all terms insecond. The specificationfirst*second indicates the cross offirst andsecond. This is the same asfirst + second + first:second.

The terms in the formula will be re-ordered so that main effects come first, followed by the interactions, all second-order, all third-order and so on: to avoid this, pass aterms object as the formula.

If you choosebinomial.glmm as thefamily.glmm, then your response should be a two-column matrix: the first column reports the number of successes and the second reports the number of failures.

The random effects for glmm are also specified symbolically. The random effects model specification is typically a list. Each element of therandom list has the formresponse ~ 0 + term. The 0 centers the random effects at 0. If you want your random effects to have a nonzero mean, then include that term in the fixed effects. Each variance component must have its own formula in the list.

To set some variance components equal to one another, use thevarcomps.equal argument. The argumentvarcomps.equal should be a vector whose length is equal to the length of the random effects list. The vector should contain positive integers, and the first element of thevarcomps.equal should be 1. To set variance components equal to one another, assign the same integer to the corresponding elements ofvarcomps.equal. For example, to set the first and second variance components equal to each other, the first two elements ofvarcomps.equal should be 1. Ifvarcomps.equal is omitted, then the variance components are assumed to be distinct.

Each distinct variance component should have a name. The length ofvarcomps.names should be equal to the number of distinct variance components. Ifvarcomps.equal is omitted, then the length ofvarcomps.names should be equal to the length ofrandom.

The package uses a relevance-weighted log density weighting scheme, similar to that described in Hu and Zidek (1997). This means that theweights argument functions in the following manner: a model built for a data set with three observations (obs. 1, obs. 2, obs. 3) and weighting scheme 1, 1, 2 gives the same results as the model built for a data set with four observations (obs. 1, obs. 2, obs. 3, obs. 3 - the last two observations are identical) and weighting scheme 1, 1, 1, 1.

Monte Carlo likelihood approximation relies on an importance sampling distribution. Though infinitely many importance sampling distributions should yield the correct MCMLEs eventually, the importance sampling distribution used in this package was chosen to reduce the computation cost. WhendoPQL isTRUE, the importance sampling distribution relies on PQL estimates (as calculated in this package). WhendoPQL isFALSE, the random effect estimates in the distribution are taken to be 0, the fixed effect estimates are taken to be 0, and the variance component estimates are taken to be 1.

This package's importance sampling distribution is a mixture of three distributions: a t centered at 0 with scale matrix determined by the PQL estimates of the variance components and withzeta degrees of freedom, a normal distribution centered at the PQL estimates of the random effects and with a variance matrix containing the PQL estimates of the variance components, and a normal distribution centered at the PQL estimates of the random effects and with a variance matrix based on the Hessian of the penalized log likelihood. The first component is included to guarantee the gradient of the MCLA has a central limit theorem. The second component is included to mirror our best guess of the distribution of the random effects. The third component is included so that the numerator and the denominator are similar when calculating the MCLA value.

The Monte Carlo sample sizem should be chosen as large as possible. You may want to run the model a couple times to begin to understand the variability inherent to Monte Carlo. There are no hard and fast rules for choosingm, and more research is needed on this area. For a general idea, I believe theBoothHobert model produces stable enough estimates atm=10^3 and thesalamander model produces stable enough estimates atm=10^5, as long asdoPQL isTRUE.

To decrease the computation time involved with larger Monte Carlo sample sizes, a parallel computing component has been added to this package. Generally, the larger the Monte Carlo sample sizem, the greater the benefit of utilizing additional cores. By default,glmm will create a cluster that uses a single core. This forces all computations to be done sequentially rather than simultaneously.

To see the summary of the model, use summary().

Value

glmm returns an object of classglmm is a list containing at least the following components:

The functionsummary (i.e.,summary.glmm) canbe used to obtain or print a summary of the results. The generic accessor functioncoef (i.e.,coef.glmm) can be used to extract the coefficients.

Author(s)

Christina Knudson

References

Geyer, C. J. (1994) On the convergence of Monte Carlo maximum likelihood calculations.Journal of the Royal Statistical Society, Series B,61, 261–274.doi:10.1111/j.2517-6161.1994.tb01976.x.

Geyer, C. J. and Thompson, E. (1992) Constrained Monte Carlo maximum likelihood for dependent data.Journal of the Royal Statistical Society, Series B,54, 657–699.doi:10.1111/j.2517-6161.1992.tb01443.x.

Knudson, C. (2016). Monte Carlo likelihood approximation for generalized linear mixed models. PhD thesis, University of Minnesota.http://hdl.handle.net/11299/178948

Knudson, C., Benson, S., Geyer, J., Jones, G. (2021) Likelihood-based inference for generalized linear mixed models: Inference with the R package glmm.Stat, 10:e339.doi:10.1002/sta4.339.

Sung, Y. J. and Geyer, C. J. (2007) Monte Carlo likelihood inference for missing data models.Annals of Statistics,35, 990–1011.doi:10.1214/009053606000001389.

Hu, F. and Zidek, J. V. (2001) The relevance weighted likelihood withapplications. In: Ahmed, S. E. and Reid N. (eds).Empirical Bayes and Likelihood Inference.Lecture Notes in Statistics,148, 211–235. Springer, New York.doi:10.1007/978-1-4613-0141-7_13.

Examples

#First, using the basic Booth and Hobert dataset#to fit a glmm with a logistic link, one variance component,#one fixed effect, and an intercept of 0. The Monte Carlo#sample size is 100 to save time.library(glmm)data(BoothHobert)set.seed(1234)mod.mcml1 <- glmm(y~0+x1, list(y~0+z1), varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=100, doPQL=TRUE)mod.mcml1$betamod.mcml1$nusummary(mod.mcml1)coef(mod.mcml1)# This next example has crossed random effects but we assume# all random effects have the same variance. The Monte Carlo#sample size is 100 to save time.#data(salamander)#set.seed(1234)#onenu <- glmm(Mate~Cross, random=list(~0+Female,~0+Male), #varcomps.names=c("only"), varcomps.equal=c(1,1), data=salamander, #family.glmm=bernoulli.glmm, m=100, debug=TRUE, doPQL=TRUE)#summary(onenu)

Monte Carlo Log Likelihood

Description

A function that calculates the Monte Carlo log likelihood evaluated at the the Monte Carlo maximum likelihood estimates returned fromglmm.

Usage

## S3 method for class 'glmm'logLik(object,...)

Arguments

Value

Author(s)

Christina Knudson

See Also

glmm for model fitting.

Examples

library(glmm)data(BoothHobert)set.seed(1234)mod <- glmm(y~0+x1, list(y~0+z1), varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=100)logLik(mod)

Monte Carlo Standard Error

Description

A function that calculates the Monte Carlo standard error for the Monte Carlo maximum likelihood estimates returned fromglmm.

Usage

mcse(object)

Arguments

Details

With maximum likelihood performed by Monte Carlo likelihood approximation, there are two sources of variability: there is variability from sample to sample and from Monte Carlo sample (of generated random effects) to Monte Carlo sample. The first source of variability (from sample to sample) is measured using standard error, which appears with the point estimates in thesummary tables. The second source of variability is due to the Monte Carlo randomness, and this is measured by the Monte Carlo standard error.

A large Monte Carlo standard error indicates the Monte Carlo sample sizem is too small.

Value

Author(s)

Christina Knudson

References

Geyer, C. J. (1994) On the convergence of Monte Carlo maximum likelihood calculations.Journal of the Royal Statistical Society, Series B,61, 261–274.

Knudson, C. (2016) Monte Carlo likelihood approximation for generalized linear mixed models. PhD thesis, University of Minnesota.http://hdl.handle.net/11299/178948

See Also

glmm for model fitting.

Examples

library(glmm)data(BoothHobert)set.seed(1234)mod <- glmm(y~0+x1, list(y~0+z1), varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=100, doPQL=TRUE)mcse(mod)

Monte Carlo Variance Covariance Matrix

Description

A function that calculates the Monte Carlo variance covariance matrix for the Monte Carlo maximum likelihood estimates returned fromglmm.

Usage

mcvcov(object)

Arguments

Details

With maximum likelihood performed by Monte Carlo likelihood approximation, there are two sources of variability: there is variability from sample to sample and from Monte Carlo sample (of generated random effects) to Monte Carlo sample. The first source of variability (from sample to sample) is measured using standard error, which appears with the point estimates in thesummary tables. The second source of variability is due to the Monte Carlo randomness, and this is measured by the Monte Carlo standard error.

A large entry in Monte Carlo variance covariance matrix indicates the Monte Carlo sample sizem is too small.

The square root of the diagonal elements represents the Monte Carlo standard errors. The off-diagonal entries represent the Monte Carlo covariance.

Value

Author(s)

Christina Knudson

References

Geyer, C. J. (1994) On the convergence of Monte Carlo maximum likelihood calculations.Journal of the Royal Statistical Society, Series B,61, 261–274.doi:10.1111/j.2517-6161.1994.tb01976.x.

Knudson, C. (2016) Monte Carlo likelihood approximation for generalized linear mixed models. PhD thesis, University of Minnesota.http://hdl.handle.net/11299/178948

See Also

glmm for model fitting.

Examples

library(glmm)data(BoothHobert)set.seed(1234)mod <- glmm(y~0+x1, list(y~0+z1), varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=100, doPQL=TRUE)mcvcov(mod)

Number of Homicide Victims Known

Description

Subjects responded to the question 'Within the past 12 months, how many people have you known personally that were victims of homicide?'

Usage

data(murder)

Format

A data frame with the following columns:

y

The number of homicide victims known personally by the subject.

race

a factor with levelsblack andwhite.

black

a dummy variable to indicate whether the subject was black.

white

a dummy variable to indicate whether the subject was white.

References

Agresti, A. (2002)Categorical Data Analysis, Second edition. Wiley.

Examples

data(murder)

Functions for the Poisson family.

Description

Given a scalareta, this calculates the cumulant and two derivatives for the Poisson family. Also checks that the data are entered correctly.

Usage

poisson.glmm()

Value

Note

This function is to be used by theglmm command.

Author(s)

Christina Knudson

See Also

glmm

Examples

poisson.glmm()$family.glmmpoisson.glmm()$cum(2)poisson.glmm()$cp(2)poisson.glmm()$cpp(2)

Radish count data set

Description

Data on life history traits for the invasive California wildradishRaphanus sativus.

Usage

data(radish2)

Format

A data frame with the following columns:

Site

Categorical. Experimental site where plant was grown. Two sites in this dataset.

Block

Categorical. Blocked nested within site.

Region

Categorical. Region from which individuals were obtained: northern, coastal California (N) or southern, inland California (S).

Pop

Categorical. Wild population nested within region.

varb

Categorical. Gives node of graphical model corresponding to each component of resp. This is useful for life history analysis (see aster package).

resp

Response vector.

id

Categorical. Indicates individual plants.

References

These data are a subset of data previously analyzed using fixed effectaster methods (R functionaster) in the following.

Ridley, C. E. and Ellstrand, N. C. (2010).Rapid evolution of morphology and adaptive life history inthe invasive California wild radish (Raphanus sativus) andthe implications for management.Evolutionary Applications,3, 64–76.

These data are a subset of data previously analyzed using random effectaster methods (R functionreaster) in the following.

Geyer, C. J., Ridley, C. E., Latta, R. G., Etterson, J. R.,and Shaw, R. G. (2013)Local Adaptation and Genetic Effects on Fitness: Calculationsfor Exponential Family Models with Random Effects.Annals of Applied Statistics,7, 1778–1795.doi:10.1214/13-AOAS653.

Examples

data(radish2)

Salamander mating data set from McCullagh and Nelder (1989)

Description

This data set presents the outcome of an experiment conducted at the University of Chicago in 1986 to study interbreeding between populations of mountain dusky salamanders (McCullagh and Nelder, 1989, Section 14.5).

Usage

data(salamander)

Format

A data frame with the following columns:

Mate

Whether the salamanders mated (1) or did not mate (0).

Cross

Cross between female and male type. A factor with four levels:R/R,R/W,W/R, andW/W. The type of the female salamander is listed first and the male is listed second. Rough Butt is represented by R and White Side is represented by W. For example,Cross=W/R indicates a White Side female was crossed with a Rough Butt male.

Male

Identification number of the male salamander. A factor.

Female

Identification number of the female salamander. A factor.

References

McCullagh P. and Nelder, J. A. (1989)Generalized Linear Models. Chapman and Hall/CRC.

Examples

data(salamander)

Standard Error

Description

A function that calculates the standard error for the Monte Carlo maximum likelihood estimates returned fromglmm.

Usage

se(object)

Arguments

Details

With maximum likelihood performed by Monte Carlo likelihood approximation, there are two sources of variability: there is variability from sample to sample and from Monte Carlo sample (of generated random effects) to Monte Carlo sample. The first source of variability (from sample to sample) is measured using standard error, which appears with the point estimates in thesummary tables. The second source of variability is due to the Monte Carlo randomness, and this is measured by the Monte Carlo standard error.

Value

Author(s)

Christina Knudson

See Also

glmm for model fitting.

mcse for calculating Monte Carlo standard error.

Examples

library(glmm)data(BoothHobert)set.seed(1234)mod <- glmm(y~0+x1, list(y~0+z1), varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=100, doPQL=TRUE)se(mod)

Summarizing GLMM Fits

Description

"summary" method for classglmm objects.

Usage

## S3 method for class 'glmm'summary(object, ...)## S3 method for class 'summary.glmm'print(x, digits = max(3, getOption("digits") - 3),      signif.stars = getOption("show.signif.stars"), ...)

Arguments

Value

The functionsummary.glmm computes and returns a list of summary statistics of the fitted generalized linear mixed model given inobject, using the components (list elements)"call" and"terms" from its argument, plus

Author(s)

Christina Knudson

See Also

The model fitting functionglmm, the genericsummary, and the functioncoefthat extracts the fixed effect coefficients.

Extract Model Variance Components

Description

A function that extracts the variance components returned fromglmm.

Usage

varcomps(object,...)

Arguments

Value

Author(s)

Christina Knudson

See Also

glmm for model fitting.coef.glmm for fixed effects coefficients.

Examples

library(glmm)data(BoothHobert)set.seed(1234)mod <- glmm(y~0+x1, list(y~0+z1), varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=100, doPQL=TRUE)varcomps(mod)

Variance-Covariance Matrix

Description

A function that calculates the variance-covariance matrix for the Monte Carlo maximum likelihood estimates returned fromglmm.

Usage

## S3 method for class 'glmm'vcov(object,...)

Arguments

Value

Author(s)

Christina Knudson

See Also

glmm for model fitting.

Examples

library(glmm)data(BoothHobert)set.seed(1234)mod <- glmm(y~0+x1, list(y~0+z1), varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=100, doPQL=TRUE)vcov(mod)