Movatterモバイル変換

Type:	Package
Title:	Linear Mixed Models with Sparse Matrix Methods and Smoothing
Description:	Provides tools for fitting linear mixed models using sparse matrix methods and variance component estimation. Applications include spline-based modeling of spatial and temporal trends using penalized splines (Boer, 2023) <doi:10.1177/1471082X231178591>.
Version:	1.0.12
Date:	2025-12-03
License:	GPL-3
Encoding:	UTF-8
LazyData:	true
Depends:	R (≥ 3.6)
Imports:	Matrix, methods, Rcpp (≥ 0.10.4), spam, splines
LinkingTo:	Rcpp
RoxygenNote:	7.3.3
Suggests:	rmarkdown, knitr, tinytest, tidyr, ggplot2, maps, sf
VignetteBuilder:	knitr
URL:	https://biometris.github.io/LMMsolver/index.html,https://github.com/Biometris/LMMsolver/
BugReports:	https://github.com/Biometris/LMMsolver/issues
NeedsCompilation:	yes
Packaged:	2025-12-05 07:49:25 UTC; rossu027
Author:	Martin Boer [aut], Bart-Jan van Rossum [aut, cre]
Maintainer:	Bart-Jan van Rossum <bart-jan.vanrossum@wur.nl>
Repository:	CRAN
Date/Publication:	2025-12-05 08:20:02 UTC

Package LMMsolver

Description

Linear Mixed Model Solver using sparse matrix algebra.

Details

An efficient and flexible system to solve sparse mixed modelequations, for models that are often used in statistical genetics.Important applications are the use of splines to model spatial or temporaltrends. Another application area is mixed model QTL analysis formultiparental populations, allowing for heterogeneous residual variance andrandom design matrices with Identity-By-Descent (IBD) probabilities.

Author(s)

Martin Boermartin.boer@wur.nl

Bart-Jan van Rossumbart-jan.vanrossum@wur.nl (maintainer)

References

Martin P. Boer (2023).Tensor product P-splines using a sparse mixed model formulation,Statistical Modelling, 23, p. 465 - 479.doi:10.1177/1471082X231178591

See Also

Useful links:

• https://biometris.github.io/LMMsolver/index.html

• https://github.com/Biometris/LMMsolver/

• Report bugs athttps://github.com/Biometris/LMMsolver/issues

construct object for Automated Differentiation Cholesky decomposition

Description

Construct object for reverse Automated Differentiation of Cholesky decomposition,with as input a list of semi-positive symmetric sparse matricesP_i, each ofdimensionq \times q. The functionADchol calculates the matrixC, the sumthe precision matricesP_i:C = \sum_{i} P_i. Next, it calculates the CholeskyDecomposition using the multiple minimum degree (MMD) algorithmof thespam package.

Usage

ADchol(lP)

Arguments

Value

An object of classADchol. This object is used to calculate the partialpartial derivatives oflog|C| in an efficient way.

References

Furrer, R., & Sain, S. R. (2010). spam: A sparse matrix R package with emphasison MCMC methods for Gaussian Markov random fields.Journal of Statistical Software, 36, 1-25.

Simulated Biomass as function of time using APSIM wheat.

Description

Simulated Biomass as function of time using APSIM wheat.

Usage

APSIMdat

Format

A data.frame with 121 rows and 4 columns.

env

Environment, Emerald in 1993

geno

Simulated genotype g001

das

Days after sowing

biomass

Simulated biomass using APSIM; medium measurement error added

References

Bustos-Korts et al. (2019) Combining Crop Growth Modeling andStatistical Genetic Modeling to Evaluate Phenotyping Strategiesdoi:10.3389/FPLS.2019.01491

Construct design matrix for B-Splines

Description

Construct design matrix for B-Splines.

Usage

Bsplines(knots, x, deriv = 0)

Arguments

Solve Linear Mixed Models

Description

Solve Linear Mixed Models using REML.

Usage

LMMsolve(  fixed,  random = NULL,  spline = NULL,  group = NULL,  ginverse = NULL,  weights = NULL,  data,  residual = NULL,  family = gaussian(),  offset = 0,  tolerance = 1e-06,  trace = FALSE,  maxit = 250,  theta = NULL,  grpTheta = NULL)

Arguments

Details

A Linear Mixed Model (LMM) has the form

y = X \beta + Z u + e, u \sim N(0,G), e \sim N(0,R)

wherey is a vector of observations,\beta is a vector with the fixedeffects,u is a vector with the random effects, ande a vector ofrandom residuals.X andZ are design matrices.

LMMsolve can fit models where the matricesG^{-1} andR^{-1} area linear combination of precision matricesQ_{G,i} andQ_{R,i}:

G^{-1} = \sum_{i} \psi_i Q_{G,i} \;, R^{-1} = \sum_{i} \phi_i Q_{R,i}

where the precision parameters\psi_i and\phi_i are estimatedusing REML. For most standard mixed models1/{\psi_i} are the variancecomponents and1/{\phi_i} the residual variances. We use a formulationin terms of precision parameters to allow for non-standard mixed models usingtensor product splines.

Value

An object of classLMMsolve representing the fitted model.SeeLMMsolveObject for a full description of the components inthis object.

See Also

LMMsolveObject,spl1D,spl2D,spl3D

Examples

## Fit models on oats.datadata(oats.data)## Fit simple model with only fixed effects.LMM1 <- LMMsolve(fixed = yield ~ rep + gen,                data = oats.data)## Fit the same model with genotype as random effect.LMM1_rand <- LMMsolve(fixed = yield ~ rep,                     random = ~gen,                     data = oats.data)## Fit the model with a 1-dimensional spline at the plot level.LMM1_spline <- LMMsolve(fixed = yield ~ rep + gen,                       spline = ~spl1D(x = plot, nseg = 20),                       data = oats.data)## Fit models on multipop data included in the package.data(multipop)## The residual variances for the two populations can be different.## Allow for heterogeneous residual variances using the residual argument.LMM2 <- LMMsolve(fixed = pheno ~ cross,                residual = ~cross,                data = multipop)## QTL-probabilities are defined by the columns pA, pB, pC.## They can be included in the random part of the model by specifying the## group argument and using grp() in the random part.# Define groups by specifying columns in data corresponding to groups in a list.# Name used in grp() should match names specified in list.lGrp <- list(QTL = 3:5)LMM2_group <- LMMsolve(fixed = pheno ~ cross,                      group = lGrp,                      random = ~grp(QTL),                      residual = ~cross,                      data = multipop)

Fitted LMMsolve Object

Description

An object of classLMMsolve returned by the LMMsolve function,representing a fitted linear mixed model. Objects of this class havemethods for the generic functions coef, fitted, residuals, loglik anddeviance.

Value

An object of classLMMsolve contains the following components:

.

Construct equally placed knots

Description

Construct equally placed knots.

Usage

PsplinesKnots(xmin, xmax, degree, nseg, cyclic = FALSE)

Arguments

Value

A numerical vector of knot positions.

Row-wise kronecker product

Description

Row-wise kronecker product

Usage

RowKronecker(X1, X2)

Arguments

Value

The row-wise kronecker product of X1 and X2.

Sea Surface Temperature

Description

Sea Surface Temperature

Usage

SeaSurfaceTemp

Format

A data.frame with 15607 rows and 4 columns.

lon

longitude

lat

latitude

sst

sea surface temperature in Kelvin

type

defines training and test set

References

Cressie et al. (2022) Basis-function models in spatial statistics.Annual Review of Statistics and Its Application.doi:10.1146/annurev-statistics-040120-020733

Uniformity trial of barley

Description

Uniformity trial of barley

Usage

barley.uniformity.trial

Format

A data.frame with 1076 rows and 3 columns

row

row coordinate

col

column coordinate

yield

yield per plot

Source

H. P. Piepho & E. R. Williams (2010). Linear variance modelsfor plant breeding trials. Plant Breeding, 129, 1-8.doi:10.1111/j.1439-0523.2009.01654.x

References

Piepho, Hans‐Peter, Martin P. Boer, and Emlyn R. Williams."Two‐dimensional P‐spline smoothing for spatial analysis of plant breeding trials."Biometrical Journal 64, no. 5 (2022): 835-857.

Standard errors for predictions

Description

Calculates the standard errors for predictionsD \hat{u},see Welham et al. 2004 and Gilmour et al. 2004 for details.

Usage

calcStandardErrors(C, D)

Arguments

Details

The prediction error variance is given byD C^{-1} D',whereC is the mixed model coefficient matrix, andD defineslinear combinations of fixed and random effects.The standard errors are given by the the square root ofthe diagonal. To calculate the standard errors in an efficient way we use that

\frac{\partial log|C + \xi_i d_i d_i'|}{\partial \xi_i} |_{\xi_i=0} = trace(C^{-1} d_i d_i') =trace(d_i' C^{-1} d_i) = d_i' C^{-1} d_i,

whered_i is rowi of matrixD. The values ofd_i' C^{-1} d_i can be calculated more efficient, avoiding thecalculation of the inverse ofC, by using Automated Differentiationof the Choleksy algorithm, see section 2.3 in Smith (1995) for details.

Value

a vector with standard errors for predictionsD \hat{u}.

References

Welham, S., Cullis, B., Gogel, B., Gilmour, A., & Thompson, R. (2004).Prediction in linear mixed models.Australian & New Zealand Journal of Statistics, 46(3), 325-347.

Smith, S. P. (1995). Differentiation of the Cholesky algorithm.Journal of Computational and Graphical Statistics, 4(2), 134-147.

Gilmour, A., Cullis, B., Welham, S., Gogel, B., & Thompson, R. (2004).An efficient computing strategy for prediction in mixed linear models.Computational statistics & data analysis, 44(4), 571-586.

Coefficients from the mixed model equations of an LMMsolve object.

Description

Obtain the coefficients from the mixed model equations of an LMMsolve object.

Usage

## S3 method for class 'LMMsolve'coef(object, se = FALSE, ...)

Arguments

Value

A list of vectors, containing the estimated effects for each fixedeffect and the predictions for each random effect in the defined linearmixed model.

Examples

## Fit model on oats datadata(oats.data)## Fit simple model with only fixed effects.LMM1 <- LMMsolve(fixed = yield ~ rep + gen,                data = oats.data)## Obtain coefficients.coefs1 <- coef(LMM1)## Obtain coefficients with standard errors.coefs2 <- coef(LMM1, se = TRUE)

Helper function for constructing Rinv

Description

Helper function for constructing Rinv

Usage

constructRinv(df, residual, weights)

Deviance of an LMMsolve object

Description

Obtain the deviance of a model fitted using LMMsolve.

Usage

## S3 method for class 'LMMsolve'deviance(object, relative = TRUE, includeConstant = TRUE, ...)

Arguments

Value

The deviance of the fitted model.

Examples

## Fit model on oats.datadata(oats.data)## Fit simple model with only fixed effects.LMM1 <- LMMsolve(fixed = yield ~ rep + gen,                data = oats.data)## Obtain deviance.deviance(LMM1)

Give diagnostics for mixed model coefficient matrix C and the choleskydecomposition

Description

Give diagnostics for mixed model coefficient matrix C and the choleskydecomposition

Usage

diagnosticsMME(object)

Arguments

Value

A summary of the mixed model coefficient matrix and its choleskidecomposition.

Examples

## Fit model on oats datadata(oats.data)## Fit simple model with only fixed effects.LMM1 <- LMMsolve(fixed = yield ~ rep + gen,                data = oats.data)## Obtain deviance.diagnosticsMME(LMM1)

Display the sparseness of the mixed model coefficient matrix

Description

Display the sparseness of the mixed model coefficient matrix

Usage

displayMME(object, cholesky = FALSE)

Arguments

Value

A plot of the sparseness of the mixed model coefficient matrix.

Examples

## Fit model on oats datadata(oats.data)## Fit simple model with only fixed effects.LMM1 <- LMMsolve(fixed = yield ~ rep + gen,                data = oats.data)## Obtain deviance.displayMME(LMM1)

Function to get the Effective Dimensions.

Description

Function to get the Effective Dimensions.

Usage

effDim(object)

Arguments

Value

A data.frame with the effective dimensions and penalties.

#' @examples## Fit model on oats datadata(oats.data)

## Fit a model with a 1-dimensional spline at the plot level.obj <- LMMsolve(fixed = yield ~ rep + gen,spline = ~spl1D(x = plot, nseg = 20),data = oats.data)effDim(obj)

Fitted values of an LMMsolve object.

Description

Obtain the fitted values from a mixed model fitted using LMMSolve.

Usage

## S3 method for class 'LMMsolve'fitted(object, ...)

Arguments

Value

A vector of fitted values.

Examples

## Fit model on oats datadata(oats.data)## Fit simple model with only fixed effects.LMM1 <- LMMsolve(fixed = yield ~ rep + gen,                data = oats.data)## Obtain fitted values.fitted1 <- fitted(LMM1)

Log-likelihood of an LMMsolve object

Description

Obtain the Restricted Maximum Log-Likelihood of a model fitted usingLMMsolve.

Usage

## S3 method for class 'LMMsolve'logLik(object, includeConstant = TRUE, ...)

Arguments

Value

The restricted maximum log-likelihood of the fitted model.

Examples

## Fit model on oats datadata(oats.data)## Fit simple model with only fixed effects.LMM1 <- LMMsolve(fixed = yield ~ rep + gen,                data = oats.data)## Obtain log-likelihood.logLik(LMM1)## Obtain log-likelihood without constant.logLik(LMM1, includeConstant = FALSE)

Function to obtain restricted log-likelihood and the first derivatives of the log-likelihood,given values for the penalty parameters

Description

Function to obtain restricted log-likelihood and the first derivatives of the log-likelihood,given values for the penalty parameters

Usage

mLogLik(object, theta)

Arguments

Value

A data.frame with logL and the first derivatives of log-likelihood

Family Object for Multinomial Model

Description

The Multinomial model is not part of the standard family. The implementationis based on Chapter 6 in Fahrmeir et al. (2013).

Usage

multinomial()

Value

An object of classfamilyLMMsolver with the following components:

References

Fahrmeir, Ludwig, Thomas Kneib, Stefan Lang, Brian Marx, Regression models.Springer Berlin Heidelberg, 2013.

Simulated QTL mapping data set

Description

Simulated QTL mapping data set

Usage

multipop

Format

A data.frame with 180 rows and 6 columns.

cross

Cross ID, two populations, AxB and AxC

ind

Genotype ID

pA

Probability that individual has alleles from parent A

pB

Probability that individual has alleles from parent B

pC

Probability that individual has alleles from parent C

pheno

Simulated phenotypic value

Alpha lattice design of spring oats

Description

Alpha lattice design of spring oats

Usage

oats.data

Format

A data.frame with 72 rows and 7 columns

plot

plot number

rep

replicate

block

incomplete block

gen

genotype

yield

dry matter yield

row

row

col

column

Details

The response is grain yield in kg per hectare. The design was an alpha designwith 24 varieties, three replicates and six incomplete blocks of size four per replicate.The 72 plots were arranged in a single linear array.

Source

J. A. John & E. R. Williams (1995). Cyclic and computergenerated designs. Chapman and Hall, London. Page 146.

References

Boer, Martin P., Hans-Peter Piepho, and Emlyn R. Williams."Linear variance, P-splines and neighbour differences for spatial adjustmentin field trials: how are they related?." JABES 25, no. 4 (2020): 676-698.

Obtain Smooth Trend.

Description

Obtain the smooth trend for models fitted with a spline component.

Usage

obtainSmoothTrend(  object,  grid = NULL,  newdata = NULL,  deriv = 0,  includeIntercept = FALSE,  which = 1)

Arguments

Value

A data.frame with predictions for the smooth trend on the specifiedgrid. The standard errors are saved if 'deriv' has default value 0.

Examples

## Fit model on oats datadata(oats.data)## Fit a model with a 1-dimensional spline at the plot level.LMM1_spline <- LMMsolve(fixed = yield ~ rep + gen,                       spline = ~spl1D(x = plot, nseg = 20),                       data = oats.data)## Obtain the smooth trend for the fitted model on a dense grid.smooth1 <- obtainSmoothTrend(LMM1_spline,                            grid = 100)## Obtain the smooth trend on a new data set - plots 10 to 40.newdat <- data.frame(plot = 10:40)smooth2 <- obtainSmoothTrend(LMM1_spline,                            newdata = newdat)## The first derivative of the smooth trend can be obtained by setting deriv = 1.smooth3 <- obtainSmoothTrend(LMM1_spline,                            grid = 100,                            deriv = 1)## For examples of higher order splines see the vignette.

Predict function

Description

Predict function

Usage

## S3 method for class 'LMMsolve'predict(  object,  newdata,  type = c("response", "link"),  se.fit = FALSE,  deriv = NULL,  ...)

Arguments

Value

A data.frame with predictions for the smooth trend on the specifiedgrid. The standard errors are saved if 'se.fit=TRUE'.

Examples

## simulate some dataf <- function(x) { 0.3 + 0.4*x + 0.2*sin(20*x) }set.seed(12)n <- 150x <- seq(0, 1, length = n)sigma2e <- 0.04y <- f(x) + rnorm(n, sd = sqrt(sigma2e))dat <- data.frame(x, y)## fit the modelobj <- LMMsolve(fixed = y ~ 1,         spline = ~spl1D(x, nseg = 50), data = dat)## make predictionsnewdat <- data.frame(x = seq(0, 1, length = 5))pred <- predict(obj, newdata = newdat, se.fit = TRUE)pred## make predictions for derivative of x:pred2 <- predict(obj, newdata = newdat, se.fit = TRUE, deriv = "x")pred2

Residuals of an LMMsolve object.

Description

Obtain the residuals from a mixed model fitted using LMMSolve.

Usage

## S3 method for class 'LMMsolve'residuals(object, ...)

Arguments

Value

A vector of residuals.

Examples

## Fit model on oats.datadata(oats.data)## Fit simple model with only fixed effects.LMM1 <- LMMsolve(fixed = yield ~ rep + gen,                data = oats.data)## Obtain fitted values.residuals1 <- residuals(LMM1)

Fit P-splines

Description

Fit multi dimensional P-splines using sparse implementation.

Usage

spl1D(  x,  nseg,  pord = 2,  degree = 3,  cyclic = FALSE,  scaleX = TRUE,  xlim = range(x),  cond = NULL,  level = NULL)spl2D(  x1,  x2,  nseg,  pord = 2,  degree = 3,  cyclic = c(FALSE, FALSE),  scaleX = TRUE,  x1lim = range(x1),  x2lim = range(x2),  cond = NULL,  level = NULL)spl3D(  x1,  x2,  x3,  nseg,  pord = 2,  degree = 3,  scaleX = TRUE,  x1lim = range(x1),  x2lim = range(x2),  x3lim = range(x3))

Arguments

Value

A list with the following elements:

• X - design matrix for fixed effect. The intercept is not included.

• Z - design matrix for random effect.

• lGinv - a list of precision matrices

• knots - a list of vectors with knot positions

• dim.f - the dimensions of the fixed effect.

• dim.r - the dimensions of the random effect.

• term.labels.f - the labels for the fixed effect terms.

• term.labels.r - the labels for the random effect terms.

• x - a list of vectors for the spline variables.

• pord - the order of the penalty.

• degree - the degree of the B-spline basis.

• scaleX - logical indicating if the fixed effects are scaled.

• EDnom - the nominal effective dimensions.

Functions

• spl2D(): 2-dimensional splines

• spl3D(): 3-dimensional splines

See Also

LMMsolve

Examples

## Fit model on oats datadata(oats.data)## Fit a model with a 1-dimensional spline at the plot level.LMM1_spline <- LMMsolve(fixed = yield ~ rep + gen,                       spline = ~spl1D(x = plot, nseg = 20),                       data = oats.data)summary(LMM1_spline)## Fit model on US precipitation data from spam package.data(USprecip, package = "spam")## Only use observed dataUSprecip <- as.data.frame(USprecip)USprecip <- USprecip[USprecip$infill == 1, ]## Fit a model with a 2-dimensional P-spline.LMM2_spline <- LMMsolve(fixed = anomaly ~ 1,                       spline = ~spl2D(x1 = lon, x2 = lat, nseg = c(41, 41)),                       data = USprecip)summary(LMM2_spline)

Summarize Linear Mixed Model fits

Description

Summary method for class "LMMsolve". Creates either a table of effectivedimensions (which = "dimensions") or a table of variances (which ="variances").

Usage

## S3 method for class 'LMMsolve'summary(object, which = c("dimensions", "variances"), ...)## S3 method for class 'summary.LMMsolve'print(x, ...)

Arguments

Value

A data.frame with either effective dimensions or variances dependingon which.

Methods (by generic)

• print(summary.LMMsolve): print summary

Examples

## Fit model on oats data.data(oats.data)## Fit simple model with only fixed effects.LMM1 <- LMMsolve(fixed = yield ~ rep + gen,                data = oats.data)## Obtain table of effective dimensions.summ1 <- summary(LMM1)print(summ1)## Obtain table of variances.summ2 <- summary(LMM1,                which = "variances")print(summ2)