Movatterモバイル変換

lava: Latent Variable Models

Description

A general implementation of Structural Equation Models with latent variables (MLE, 2SLS, and composite likelihood estimators) with both continuous, censored, and ordinal outcomes (Holst and Budtz-Joergensen (2013)doi:10.1007/s00180-012-0344-y). Mixture latent variable models and non-linear latent variable models (Holst and Budtz-Joergensen (2020)doi:10.1093/biostatistics/kxy082). The package also provides methods for graph exploration (d-separation, back-door criterion), simulation of general non-linear latent variable models, and estimation of influence functions for a broad range of statistical models.

A general implementation of Structural Equation Models wth latent variables(MLE, 2SLS, and composite likelihood estimators) with both continuous,censored, and ordinal outcomes (Holst and Budtz-Joergensen (2013)<doi:10.1007/s00180-012-0344-y>). Mixture latent variable models andnon-linear latent variable models (Holst and Budtz-Joergensen (2020)<doi:10.1093/biostatistics/kxy082>). The package also provides methods forgraph exploration (d-separation, back-door criterion), simulation of generalnon-linear latent variable models, and estimation of influence functions fora broad range of statistical models.

Author(s)

Maintainer: Klaus K. Holstklaus@holst.it

Other contributors:

• Brice Ozenne [contributor]

• Thomas Gerds [contributor]

Klaus K. Holst Maintainer: <klaus@holst.it>

See Also

Useful links:

• https://kkholst.github.io/lava/

• Report bugs athttps://github.com/kkholst/lava/issues

Examples

lava()

Concatenation operator

Description

For matrices a block-diagonal matrix is created. For all otherdata types he operator is a wrapper ofpaste.

Usage

x %++% y

Arguments

Details

Concatenation operator

Author(s)

Klaus K. Holst

See Also

blockdiag,paste,cat,

Examples

## Block diagonalmatrix(rnorm(25),5)%++%matrix(rnorm(25),5)## String concatenation"Hello "%++%" World"## Function compositionf <- log %++% expf(2)

Matching operator (x not in y) oposed to the%in%-operator (x in y)

Description

Matching operator

Usage

x %ni% y

Arguments

Value

A logical vector.

Author(s)

Klaus K. Holst

See Also

match

Examples

1:10 %ni% c(1,5,10)

Apply a Function to a Data Frame Split by Factors

Description

Apply a Function to a Data Frame Split by Factors

Usage

By(x, INDICES, FUN, COLUMNS, array = FALSE, ...)

Arguments

Details

Simple wrapper of the 'by' function

Author(s)

Klaus K. Holst

Examples

By(datasets::CO2,~Treatment+Type,colMeans,~conc)By(datasets::CO2,~Treatment+Type,colMeans,~conc+uptake)

Generate a transparent RGB color

Description

This function transforms a standard color (e.g. "red") into antransparent RGB-color (i.e. alpha-blend<1).

Usage

Col(col, alpha = 0.2, locate = 0)

Arguments

Details

This only works for certain graphics devices (Cairo-X11 (x11 as of R>=2.7), quartz, pdf, ...).

Value

A character vector with elements of 7 or 9 characters, '#'followed by the red, blue, green and optionally alpha values inhexadecimal (after rescaling to '0 ... 255').

Author(s)

Klaus K. Holst

Examples

plot(runif(1000),cex=runif(1000,0,4),col=Col(c("darkblue","orange"),0.5),pch=16)

Report estimates across different models

Description

Report estimates across different models

Usage

Combine(x, ...)

Arguments

Author(s)

Klaus K. Holst

Examples

data(serotonin)m1 <- lm(cau ~ age*gene1 + age*gene2,data=serotonin)m2 <- lm(cau ~ age + gene1,data=serotonin)m3 <- lm(cau ~ age*gene2,data=serotonin)Combine(list(A=m1,B=m2,C=m3),fun=function(x)     c("_____"="",R2=" "%++%format(summary(x)$r.squared,digits=2)))

Create a Data Frame from All Combinations of Factors

Description

Create a Data Frame from All Combinations of Factors

Usage

Expand(`_data`, ...)

Arguments

Details

Simple wrapper of the 'expand.grid' function. If x is a tablethen a data frame is returned with one row pr individualobservation.

Author(s)

Klaus K. Holst

Examples

dd <- Expand(iris, Sepal.Length=2:8, Species=c("virginica","setosa"))summary(dd)T <- with(warpbreaks, table(wool, tension))Expand(T)

Extract graph

Description

Extract or replace graph object

Usage

Graph(x, ...)Graph(x, ...) <- value

Arguments

Author(s)

Klaus K. Holst

See Also

Model

Examples

m <- lvm(y~x)Graph(m)

Finds elements in vector or column-names in data.frame/matrix

Description

Pattern matching in a vector or column names of a data.frame or matrix.

Usage

Grep(x, pattern, subset = TRUE, ignore.case = TRUE, ...)

Arguments

Value

A data.frame with 2 columns with the indices in the first and thematching names in the second.

Author(s)

Klaus K. Holst

See Also

grep, andagrep for approximate stringmatching,

Examples

data(iris)head(Grep(iris,"(len)|(sp)"))

Extract i.i.d. decomposition (influence function) from model object

Description

Extract i.i.d. decomposition (influence function) from model object

Usage

IC(x,...)## Default S3 method:IC(x, bread, id=NULL, folds=0, maxsize=(folds>0)*1e6,...)

Arguments

Examples

m <- lvm(y~x+z)distribution(m, ~y+z) <- binomial.lvm("logit")d <- sim(m,1e3)g <- glm(y~x+z,data=d,family=binomial)var_ic(IC(g))

Missing value generator

Description

Missing value generator

Usage

Missing(object, formula, Rformula, missing.name, suffix = "0", ...)

Arguments

Details

This function adds a binary variable to a givenlvm modeland also a variable which is equal to the original variable wherethe binary variable is equal to zero

Value

lvm object

Author(s)

Thomas A. Gerds <tag@biostat.ku.dk>

Examples

library(lava)set.seed(17)m <- lvm(y0~x01+x02+x03)m <- Missing(m,formula=x1~x01,Rformula=R1~0.3*x02+-0.7*x01,p=0.4)sim(m,10)m <- lvm(y~1)m <- Missing(m,"y","r")## same as## m <- Missing(m,y~1,r~1)sim(m,10)## same asm <- lvm(y~1)Missing(m,"y") <- r~xsim(m,10)m <- lvm(y~1)m <- Missing(m,"y","r",suffix=".")## same as## m <- Missing(m,"y","r",missing.name="y.")## same as## m <- Missing(m,y.~y,"r")sim(m,10)

Extract model

Description

Extract or replace model object

Usage

Model(x, ...)Model(x, ...) <- value

Arguments

Value

Returns a model object (e.g.lvm ormultigroup)

Author(s)

Klaus K. Holst

See Also

Graph

Examples

m <- lvm(y~x)e <- estimate(m, sim(m,100))Model(e)

Convert to/from NA

Description

Convert vector to/from NA

Usage

NA2x(s, x = 0)

Arguments

Value

A vector with same dimension and class ass.

Author(s)

Klaus K. Holst

Examples

##' x2NA(1:10, 1:5)NA2x(x2NA(c(1:10),5),5)##'

Newton-Raphson method

Description

Newton-Raphson method

Usage

NR(  start,  objective = NULL,  gradient = NULL,  hessian = NULL,  control,  args = NULL,  ...)

Arguments

Details

control should be a list with one or more of the following components:

• trace integer for which output is printed each 'trace'th iteration

• iter.max number of iterations

• stepsize: Step size (default 1)

• nstepsize: Increase stepsize every nstepsize iteration (from stepsize to 1)

• tol: Convergence criterion (gradient)

• epsilon: threshold used in pseudo-inverse

• backtrack: In each iteration reduce stepsize unless solution is improved according to criterion (gradient, armijo, curvature, wolfe)

Examples

# Objective function with gradient and hessian as attributesf <- function(z) {   x <- z[1]; y <- z[2]   val <- x^2 + x*y^2 + x + y   structure(val, gradient=c(2*x+y^2+1, 2*y*x+1),             hessian=rbind(c(2,2*y),c(2*y,2*x)))}NR(c(0,0),f)# Parsing arguments to the function andg <- function(x,y) (x*y+1)^2NR(0, gradient=g, args=list(y=2), control=list(trace=1,tol=1e-20))

Dose response calculation for binomial regression models

Description

Dose response calculation for binomial regression models

Usage

PD(  model,  intercept = 1,  slope = 2,  prob = NULL,  x,  level = 0.5,  ci.level = 0.95,  vcov,  family,  EB = NULL)

Arguments

Author(s)

Klaus K. Holst

Generic print method

Description

Nicer print method for tabular data. Falls back to standard print method forall other data types.

Usage

Print(x, n = 5, digits = max(3, getOption("digits") - 3), ...)

Arguments

Define range constraints of parameters

Description

Define range constraints of parameters

Usage

Range.lvm(a = 0, b = 1)

Arguments

Value

function

Author(s)

Klaus K. Holst

Add variable to (model) object

Description

Generic method for adding variables to model object

Usage

addvar(x, ...)

Arguments

Author(s)

Klaus K. Holst

Backdoor criterion

Description

Check backdoor criterion of a lvm object

Usage

backdoor(object, f, cond, ..., return.graph = FALSE)

Arguments

Examples

m <- lvm(y~c2,c2~c1,x~c1,m1~x,y~m1, v1~c3, x~c3,v1~y,         x~z1, z2~z1, z2~z3, y~z3+z2+g1+g2+g3)ll <- backdoor(m, y~x)backdoor(m, y~x|c1+z1+g1)

Label elements of object

Description

Generic method for labeling elements of an object

Usage

baptize(x, ...)

Arguments

Author(s)

Klaus K. Holst

Define constant risk difference or relative risk association for binary exposure

Description

Set up model as defined in Richardson, Robins and Wang (2017).

Usage

binomial.rd(  x,  response,  exposure,  target.model,  nuisance.model,  exposure.model = binomial.lvm(),  ...)

Arguments

Combine matrices to block diagonal structure

Description

Combine matrices to block diagonal structure

Usage

blockdiag(x, ..., pad = 0)

Arguments

Author(s)

Klaus K. Holst

Examples

A <- diag(3)+1blockdiag(A,A,A,pad=NA)

Longitudinal Bone Mineral Density Data (Wide format)

Description

Bone Mineral Density Data consisting of 112 girls randomized to receivecalcium og placebo. Longitudinal measurements of bone mineral density(g/cm^2) measured approximately every 6th month in 3 years.

Format

data.frame

Source

Vonesh & Chinchilli (1997), Table 5.4.1 on page 228.

See Also

calcium

Data

Description

Description

Format

data.frame

Generic bootstrap method

Description

Generic method for calculating bootstrap statistics

Usage

bootstrap(x, ...)

Arguments

Author(s)

Klaus K. Holst

See Also

bootstrap.lvmbootstrap.lvmfit

Calculate bootstrap estimates of a lvm object

Description

Draws non-parametric bootstrap samples

Usage

## S3 method for class 'lvm'bootstrap(x,R=100,data,fun=NULL,control=list(),                          p, parametric=FALSE, bollenstine=FALSE,                          constraints=TRUE,sd=FALSE, mc.cores,                          future.args=list(future.seed=TRUE),                          ...)## S3 method for class 'lvmfit'bootstrap(x,R=100,data=model.frame(x),                             control=list(start=coef(x)),                             p=coef(x), parametric=FALSE, bollenstine=FALSE,                             estimator=x$estimator,weights=Weights(x),...)

Arguments

Value

Abootstrap.lvm object.

Author(s)

Klaus K. Holst

See Also

confint.lvmfit

Examples

m <- lvm(y~x)d <- sim(m,100)e <- estimate(lvm(y~x), data=d) ## Reduce Ex.TimingsB <- bootstrap(e,R=50,mc.cores=1)B

Simulated data

Description

Simulated data

Format

data.frame

Source

Simulated

Longitudinal Bone Mineral Density Data

Description

Bone Mineral Density Data consisting of 112 girls randomized to receivecalcium og placebo. Longitudinal measurements of bone mineral density(g/cm^2) measured approximately every 6th month in 3 years.

Format

A data.frame containing 560 (incomplete) observations. The 'person'column defines the individual girls of the study with measurements atvisiting times 'visit', and age in years 'age' at the time of visit. Thebone mineral density variable is 'bmd' (g/cm^2).

Source

Vonesh & Chinchilli (1997), Table 5.4.1 on page 228.

Generic cancel method

Description

Generic cancel method

Usage

cancel(x, ...)

Arguments

Author(s)

Klaus K. Holst

Extract children or parent elements of object

Description

Generic method for memberships from object (e.g. a graph)

Usage

children(object, ...)

Arguments

Author(s)

Klaus K. Holst

Identify points on plot

Description

Extension of theidentify function

Usage

## Default S3 method:click(x, y=NULL, label=TRUE, n=length(x), pch=19, col="orange", cex=3, ...)idplot(x, y ,..., id=list(), return.data=FALSE)

Arguments

Details

For the usual 'X11' device the identification process isterminated by pressing any mouse button other than the first. Forthe 'quartz' device the process is terminated by pressing eitherthe pop-up menu equivalent (usually second mouse button or'Ctrl'-click) or the 'ESC' key.

Author(s)

Klaus K. Holst

See Also

idplot,identify

Examples

if (interactive()) {    n <- 10; x <- seq(n); y <- runif(n)    plot(y ~ x); click(x,y)    data(iris)    l <- lm(Sepal.Length ~ Sepal.Width*Species,iris)    res <- plotConf(l,var2="Species")## ylim=c(6,8), xlim=c(2.5,3.3))    with(res, click(x,y))    with(iris, idplot(Sepal.Length,Petal.Length))}

Closed testing procedure

Description

Given p hypotheses H1, ..., Hp all 2^p-1 intersection hypotheses arecalculated and adjusted p-values are obtained for Hj is calculated as themax p-value of all intersection hypotheses containing Hj. Example, for p=3,the adjusted p-value for H1 will be obtained from {(H1, H2, H3), (H1,H2),(H1,H3), (H1)}.

Usage

closed_testing(object, test = test_wald, ...)

Arguments

Details

The function 'test' should be a function 'function(object, index,...)' which as its first argument takes an 'estimate' object and and witan argument 'index' which is a integer vector specifying whichsubcomponents of 'object' to test. The ellipsis argument can be any otherarguments used in the test function. The functiontest_wald is anexample of valid test function (which has an additional argument 'null' inreference to the above mentioned ellipsis arguments).

References

Marcus, R; Peritz, E; Gabriel, KR (1976)."On closed testing procedures with special reference to ordered analysisof variance". Biometrika. 63 (3): 655–660.

Examples

m <- lvm()regression(m, c(y1,y2,y3,y4)~x) <- c(0, 0.25, 0, 0.25)regression(m, to=endogenous(m), from="u") <- 1variance(m,endogenous(m)) <- 1set.seed(1)d <- sim(m, 200)l1 <- lm(y1~x,d)l2 <- lm(y2~x,d)l3 <- lm(y3~x,d)l4 <- lm(y4~x,d)(a <- merge(l1, l2, l3, l4, subset=2))if (requireNamespace("mets",quietly=TRUE)) {   alpha_zmax(a)}adj <- closed_testing(a)adjadj$p.valuesummary(adj)

Add color-bar to plot

Description

Add color-bar to plot

Usage

colorbar(  clut = Col(rev(rainbow(11, start = 0, end = 0.69)), alpha),  x.range = c(-0.5, 0.5),  y.range = c(-0.1, 0.1),  values = seq(clut),  digits = 2,  label.offset,  srt = 45,  cex = 0.5,  border = NA,  alpha = 0.5,  position = 1,  direction = c("horizontal", "vertical"),  ...)

Arguments

Examples

## Not run: plotNeuro(x,roi=R,mm=-18,range=5)colorbar(clut=Col(rev(rainbow(11,start=0,end=0.69)),0.5),         x=c(-40,40),y.range=c(84,90),values=c(-5:5))colorbar(clut=Col(rev(rainbow(11,start=0,end=0.69)),0.5),         x=c(-10,10),y.range=c(-100,50),values=c(-5:5),         direction="vertical",border=1)## End(Not run)

Finds the unique commutation matrix

Description

Finds the unique commutation matrix K:K vec(A) = vec(A^t)

Usage

commutation(m, n = m)

Arguments

Author(s)

Klaus K. Holst

Statistical tests

Description

Performs Likelihood-ratio, Wald and score tests

Usage

compare(object, ...)

Arguments

Value

Matrix of test-statistics and p-values

Author(s)

Klaus K. Holst

See Also

modelsearch,equivalence

Examples

m <- lvm();regression(m) <- c(y1,y2,y3) ~ eta; latent(m) <- ~etaregression(m) <- eta ~ xm2 <- regression(m, c(y3,eta) ~ x)set.seed(1)d <- sim(m,1000)e <- estimate(m,d)e2 <- estimate(m2,d)compare(e)compare(e,e2) ## LRT, H0: y3<-x=0compare(e,scoretest=y3~x) ## Score-test, H0: y3~x=0compare(e2,par=c("y3~x")) ## Wald-test, H0: y3~x=0B <- diag(2); colnames(B) <- c("y2~eta","y3~eta")compare(e2,contrast=B,null=c(1,1))B <- rep(0,length(coef(e2))); B[1:3] <- 1compare(e2,contrast=B)compare(e,scoretest=list(y3~x,y2~x))

Composite Likelihood for probit latent variable models

Description

Estimate parameters in a probit latent variable model via a compositelikelihood decomposition.

Usage

complik(  x,  data,  k = 2,  type = c("all", "nearest"),  pairlist,  messages = 0,  estimator = "normal",  quick = FALSE,  ...)

Arguments

Value

An object of classestimate.complik inheriting methods fromlvm

Author(s)

Klaus K. Holst

See Also

estimate

Examples

m <- lvm(c(y1,y2,y3)~b*x+1*u[0],latent=~u)ordinal(m,K=2) <- ~y1+y2+y3d <- sim(m,50,seed=1)if (requireNamespace("mets", quietly=TRUE)) {   e1 <- complik(m,d,control=list(trace=1),type="all")}

Add Confidence limits bar to plot

Description

Add Confidence limits bar to plot

Usage

confband(  x,  lower,  upper,  center = NULL,  line = TRUE,  delta = 0.07,  centermark = 0.03,  pch,  blank = TRUE,  vert = TRUE,  polygon = FALSE,  step = FALSE,  ...)

Arguments

Author(s)

Klaus K. Holst

See Also

confband

Examples

plot(0,0,type="n",xlab="",ylab="")confband(0.5,-0.5,0.5,0,col="darkblue")confband(0.8,-0.5,0.5,0,col="darkred",vert=FALSE,pch=1,cex=1.5)set.seed(1)K <- 20est <- rnorm(K)se <- runif(K,0.2,0.4)x <- cbind(est,est-2*se,est+2*se,runif(K,0.5,2))x[c(3:4,10:12),] <- NArownames(x) <- unlist(lapply(letters[seq(K)],function(x) paste(rep(x,4),collapse="")))rownames(x)[which(is.na(est))] <- ""signif <- sign(x[,2])==sign(x[,3])forestplot(x,text.right=FALSE)forestplot(x[,-4],sep=c(2,15),col=signif+1,box1=TRUE,delta=0.2,pch=16,cex=1.5)forestplot(x,vert=TRUE,text=FALSE)forestplot(x,vert=TRUE,text=FALSE,pch=NA)##forestplot(x,vert=TRUE,text.vert=FALSE)##forestplot(val,vert=TRUE,add=TRUE)z <- seq(10)zu <- c(z[-1],10)plot(z,type="n")confband(z,zu,rep(0,length(z)),col=Col("darkblue"),polygon=TRUE,step=TRUE)confband(z,zu,zu-2,col=Col("darkred"),polygon=TRUE,step=TRUE)z <- seq(0,1,length.out=100)plot(z,z,type="n")confband(z,z,z^2,polygon="TRUE",col=Col("darkblue"))set.seed(1)k <- 10x <- seq(k)est <- rnorm(k)sd <- runif(k)val <- cbind(x,est,est-sd,est+sd)par(mfrow=c(1,2))plot(0,type="n",xlim=c(0,k+1),ylim=range(val[,-1]),axes=FALSE,xlab="",ylab="")axis(2)confband(val[,1],val[,3],val[,4],val[,2],pch=16,cex=2)plot(0,type="n",ylim=c(0,k+1),xlim=range(val[,-1]),axes=FALSE,xlab="",ylab="")axis(1)confband(val[,1],val[,3],val[,4],val[,2],pch=16,cex=2,vert=FALSE)x <- seq(0, 3, length.out=20)y <- cos(x)yl <- y - 1yu <- y + 1plot_region(x, y, yl, yu)plot_region(x, y, yl, yu, type='s', col="darkblue", add=TRUE)

Calculate confidence limits for parameters

Description

Calculate Wald og Likelihood based (profile likelihood) confidence intervals

Usage

## S3 method for class 'lvmfit'confint(  object,  parm = seq_len(length(coef(object))),  level = 0.95,  profile = FALSE,  curve = FALSE,  n = 20,  interval = NULL,  lower = TRUE,  upper = TRUE,  ...)

Arguments

Details

Calculates either Wald confidence limits:

\hat{\theta} \pmz_{\alpha/2}*\hat\sigma_{\hat\theta}

or profile likelihood confidencelimits, defined as the set of value\tau:

logLik(\hat\theta_{\tau},\tau)-logLik(\hat\theta)< q_{\alpha}/2

whereq_{\alpha} is the\alpha fractile of the\chi^2_1distribution, and\hat\theta_{\tau} are obtained by maximizing thelog-likelihood with tau being fixed.

Value

A 2xp matrix with columns of lower and upper confidence limits

Author(s)

Klaus K. Holst

See Also

bootstrap{lvm}

Examples

m <- lvm(y~x)d <- sim(m,100)e <- estimate(lvm(y~x), d)confint(e,3,profile=TRUE)confint(e,3) ## Reduce Ex.timingsB <- bootstrap(e,R=50)B

Conformal prediction

Description

Conformal predicions using locally weighted conformal inference with a split-conformal algorithm

Usage

confpred(object, data, newdata = data, alpha = 0.05, mad, ...)

Arguments

Value

data.frame with fitted (fit), lower (lwr) and upper (upr) predictions bands.

Examples

set.seed(123)n <- 200x <- seq(0,6,length.out=n)delta <- 3ss <- exp(-1+1.5*cos((x-delta)))ee <- rnorm(n,sd=ss)y <- (x-delta)+3*cos(x+4.5-delta)+eed <- data.frame(y=y,x=x)newd <- data.frame(x=seq(0,6,length.out=50))cc <- confpred(lm(y~splines::ns(x,knots=c(1,3,5)),data=d), data=d, newdata=newd)if (interactive()) {plot(y~x,pch=16,col=lava::Col("black"),ylim=c(-10,10),xlab="X",ylab="Y")with(cc,     lava::confband(newd$x,lwr,upr,fit,        lwd=3,polygon=TRUE,col=Col("blue"),border=FALSE))}

Add non-linear constraints to latent variable model

Description

Add non-linear constraints to latent variable model

Usage

## Default S3 replacement method:constrain(x,par,args,endogenous=TRUE,...) <- value## S3 replacement method for class 'multigroup'constrain(x,par,k=1,...) <- valueconstraints(object,data=model.frame(object),vcov=object$vcov,level=0.95,                        p=pars.default(object),k,idx,...)

Arguments

Details

Add non-linear parameter constraints as well as non-linear associationsbetween covariates and latent or observed variables in the model (non-linearregression).

As an example we will specify the follow multiple regression model:

E(Y|X_1,X_2) = \alpha + \beta_1 X_1 + \beta_2 X_2

V(Y|X_1,X_2)= v

which is defined (with the appropiate parameter labels) as

m <- lvm(y ~ f(x,beta1) + f(x,beta2))

intercept(m) <- y ~ f(alpha)

covariance(m) <- y ~ f(v)

The somewhat strained parameter constraint

v =\frac{(beta1-beta2)^2}{alpha}

can then specified as

constrain(m,v ~ beta1 + beta2 + alpha) <- function(x)(x[1]-x[2])^2/x[3]

A subset of the argumentsargs can be covariates in the model,allowing the specification of non-linear regression models. As an examplethe non-linear regression model

E(Y\mid X) = \nu + \Phi(\alpha +\beta X)

where\Phi denotes the standard normal cumulativedistribution function, can be defined as

m <- lvm(y ~ f(x,0)) # No linear effect of x

Next we add three new parameters using theparameter assigmentfunction:

parameter(m) <- ~nu+alpha+beta

The intercept ofY is defined asmu

intercept(m) <- y ~ f(mu)

And finally the newly added intercept parametermu is defined as theappropiate non-linear function of\alpha,\nu and\beta:

constrain(m, mu ~ x + alpha + nu) <- function(x)pnorm(x[1]*x[2])+x[3]

Theconstraints function can be used to show the estimated non-linearparameter constraints of an estimated model object (lvmfit ormultigroupfit). Callingconstrain with no additional argumentsbeyoundx will return a list of the functions and parameter namesdefining the non-linear restrictions.

The gradient function can optionally be added as an attributegrad tothe return value of the function defined byvalue. In this case theanalytical derivatives will be calculated via the chain rule when evaluatingthe corresponding score function of the log-likelihood. If the gradientattribute is omitted the chain rule will be applied on a numericapproximation of the gradient.

Value

Alvm object.

Author(s)

Klaus K. Holst

See Also

regression,intercept,covariance

Examples

################################# Non-linear parameter constraints 1##############################m <- lvm(y ~ f(x1,gamma)+f(x2,beta))covariance(m) <- y ~ f(v)d <- sim(m,100)m1 <- m; constrain(m1,beta ~ v) <- function(x) x^2## Define slope of x2 to be the square of the residual variance of y## Estimate both restricted and unrestricted modele <- estimate(m,d,control=list(method="NR"))e1 <- estimate(m1,d)p1 <- coef(e1)p1 <- c(p1[1:2],p1[3]^2,p1[3])## Likelihood of unrestricted model evaluated in MLE of restricted modellogLik(e,p1)## Likelihood of restricted model (MLE)logLik(e1)################################# Non-linear regression################################ Simulate datam <- lvm(c(y1,y2)~f(x,0)+f(eta,1))latent(m) <- ~etacovariance(m,~y1+y2) <- "v"intercept(m,~y1+y2) <- "mu"covariance(m,~eta) <- "zeta"intercept(m,~eta) <- 0set.seed(1)d <- sim(m,100,p=c(v=0.01,zeta=0.01))[,manifest(m)]d <- transform(d,               y1=y1+2*pnorm(2*x),               y2=y2+2*pnorm(2*x))## Specify model and estimate parametersconstrain(m, mu ~ x + alpha + nu + gamma) <- function(x) x[4]*pnorm(x[3]+x[1]*x[2]) ## Reduce Ex.Timingse <- estimate(m,d,control=list(trace=1,constrain=TRUE))constraints(e,data=d)## Plot model-fitplot(y1~x,d,pch=16); points(y2~x,d,pch=16,col="gray")x0 <- seq(-4,4,length.out=100)lines(x0,coef(e)["nu"] + coef(e)["gamma"]*pnorm(coef(e)["alpha"]*x0))################################# Multigroup model################################# Define two modelsm1 <- lvm(y ~ f(x,beta)+f(z,beta2))m2 <- lvm(y ~ f(x,psi) + z)### And simulate data from themd1 <- sim(m1,500)d2 <- sim(m2,500)### Add 'non'-linear parameter constraintconstrain(m2,psi ~ beta2) <- function(x) x## Add parameter beta2 to model 2, now beta2 exists in both modelsparameter(m2) <- ~ beta2ee <- estimate(list(m1,m2),list(d1,d2),control=list(method="NR"))summary(ee)m3 <- lvm(y ~ f(x,beta)+f(z,beta2))m4 <- lvm(y ~ f(x,beta2) + z)e2 <- estimate(list(m3,m4),list(d1,d2),control=list(method="NR"))e2

Create contrast matrix

Description

Create contrast matrix typically for use with 'estimate' (Wald tests).

Usage

contr(p, n, diff = TRUE, ...)

Arguments

Examples

contr(2,n=5)contr(as.list(2:4),n=5)contr(list(1,2,4),n=5)contr(c(2,3,4),n=5)contr(list(c(1,3),c(2,4)),n=5)contr(list(c(1,3),c(2,4),5))parsedesign(c("aa","b","c"),"?","?",diff=c(FALSE,TRUE))## All pairs comparisons:pdiff <- function(n) lava::contr(lapply(seq(n-1), function(x) seq(x, n)))pdiff(4)

Generic method for extracting correlation coefficients of model object

Description

Generic correlation method

Usage

correlation(x, ...)

Arguments

Author(s)

Klaus K. Holst

Add covariance structure to Latent Variable Model

Description

Define covariances between residual terms in alvm-object.

Usage

## S3 replacement method for class 'lvm'covariance(object, var1=NULL, var2=NULL, constrain=FALSE, pairwise=FALSE,...) <- value

Arguments

Details

Thecovariance function is used to specify correlation structurebetween residual terms of a latent variable model, using a formula syntax.

For instance, a multivariate model with three response variables,

Y_1 = \mu_1 + \epsilon_1

Y_2 = \mu_2 + \epsilon_2

Y_3 = \mu_3 + \epsilon_3

can be specified as

m <- lvm(~y1+y2+y3)

Pr. default the two variables are assumed to be independent. To add acovariance parameterr = cov(\epsilon_1,\epsilon_2), we execute thefollowing code

covariance(m) <- y1 ~ f(y2,r)

The special functionf and its second argument could be omitted thusassigning an unique parameter the covariance betweeny1 andy2.

Similarily the marginal variance of the two response variables can be fixedto be identical (var(Y_i)=v) via

covariance(m) <- c(y1,y2,y3) ~ f(v)

To specify a completely unstructured covariance structure, we can call

covariance(m) <- ~y1+y2+y3

All the parameter values of the linear constraints can be given as the righthandside expression of the assigment functioncovariance<- if thefirst (and possibly second) argument is defined as well. E.g:

covariance(m,y1~y1+y2) <- list("a1","b1")

covariance(m,~y2+y3) <- list("a2",2)

Defines

var(\epsilon_1) = a1

var(\epsilon_2) = a2

var(\epsilon_3) = 2

cov(\epsilon_1,\epsilon_2) = b1

Parameter constraints can be cleared by fixing the relevant parameters toNA (see also theregression method).

The functioncovariance (called without additional arguments) can beused to inspect the covariance constraints of alvm-object.

Value

Alvm-object

Author(s)

Klaus K. Holst

See Also

regression<-,intercept<-,constrain<-parameter<-,latent<-,cancel<-,kill<-

Examples

m <- lvm()### Define covariance between residuals terms of y1 and y2covariance(m) <- y1~y2covariance(m) <- c(y1,y2)~f(v) ## Same marginal variancecovariance(m) ## Examine covariance structure

Split data into folds

Description

Split data into folds

Usage

csplit(x, p = NULL, replace = FALSE, return.index = FALSE, k = 2, ...)

Arguments

Author(s)

Klaus K. Holst

Examples

foldr(5,2,rep=2)csplit(10,3)csplit(iris[1:10,]) ## Split in two sets 1:(n/2) and (n/2+1):ncsplit(iris[1:10,],0.5)

Adds curly brackets to plot

Description

Adds curly brackets to plot

Usage

curly(  x,  y,  len = 1,  theta = 0,  wid,  shape = 1,  col = 1,  lwd = 1,  lty = 1,  grid = FALSE,  npoints = 50,  text = NULL,  offset = c(0.05, 0))

Arguments

Examples

if (interactive()) {plot(0,0,type="n",axes=FALSE,xlab="",ylab="")curly(x=c(1,0),y=c(0,1),lwd=2,text="a")curly(x=c(1,0),y=c(0,1),lwd=2,text="b",theta=pi)curly(x=-0.5,y=0,shape=1,theta=pi,text="c")curly(x=0,y=0,shape=1,theta=0,text="d")curly(x=0.5,y=0,len=0.2,theta=pi/2,col="blue",lty=2)curly(x=0.5,y=-0.5,len=0.2,theta=-pi/2,col="red",shape=1e3,text="e")}

50 patients from Monash Medical Centre, Melbourne

Description

Diagnosis of depression (DSM-III-R MDD, Dysthymia, Adjustment Disorder withDepressed Mood, Depression NOS), Beck Depression Inventory (BDI) (Beck etal., 1961) General Health Questionnaire (GHQ) (Goldberg & Williams, 1988)

Format

data.frame

Source

Clarke, D. M., Smith, G. C., & Herrman, H. E. (1993). A comparativestudy of screening instruments for mental disorders in general hospitalpatients. International Journal Psychiatry in Medicine, 23, pp. 323-337.

McKenzie et al. (1996). Comparing correlated Kappas by resampling: Is onelevel of agreement significantly different from another? J. Psychiat. Res.30 (6), pp. 483-492.

Returns device-coordinates and plot-region

Description

Returns device-coordinates and plot-region

Usage

devcoords()

Value

Alist with elements

Author(s)

Klaus K. Holst

Calculate diagnostic tests for 2x2 table

Description

Calculate prevalence, sensitivity, specificity, and positive andnegative predictive values

Usage

diagtest(  table,  positive = 2,  exact = FALSE,  p0 = NA,  confint = c("logit", "arcsin", "pseudoscore", "exact"),  ...)

Arguments

Details

Table should be in the format with outcome in columns andtest in rows. Data.frame should be with test in the firstcolumn and outcome in the second column.

Author(s)

Klaus Holst

Examples

M <- as.table(matrix(c(42,12,                       35,28),ncol=2,byrow=TRUE,                     dimnames=list(rater=c("no","yes"),gold=c("no","yes"))))diagtest(M,exact=TRUE)

Check d-separation criterion

Description

Check for conditional independence (d-separation)

Usage

## S3 method for class 'lvm'dsep(object, x, cond = NULL, return.graph = FALSE, ...)

Arguments

Details

The argument 'x' can be given as a formula, e.g. x~y|z+vor ~x+y|z+v With everything on the rhs of the bar defining thevariables on which to condition on.

Examples

m <- lvm(x5 ~ x4+x3, x4~x3+x1, x3~x2, x2~x1)if (interactive()) {plot(m,layoutType='neato')}dsep(m,x5~x1|x2+x4)dsep(m,x5~x1|x3+x4)dsep(m,~x1+x2+x3|x4)

Identify candidates of equivalent models

Description

Identifies candidates of equivalent models

Usage

equivalence(x, rel, tol = 0.001, k = 1, omitrel = TRUE, ...)

Arguments

Author(s)

Klaus K. Holst

See Also

compare,modelsearch

Estimate parameters and influence function.

Description

Estimate parameters for the sample mean, variance, and quantiles

Usage

## S3 method for class 'array'estimate(x, type = "mean", probs = 0.5, ...)

Arguments

Estimation of functional of parameters

Description

Estimation of functional of parameters. Wald tests, robust standard errors,cluster robust standard errors, LRT (whenf is not a function)...

Usage

## Default S3 method:estimate(  x = NULL,  f = NULL,  ...,  data,  id,  iddata,  stack = TRUE,  average = FALSE,  subset,  score.deriv,  level = 0.95,  IC = robust,  type = c("robust", "df", "mbn"),  keep,  use,  regex = FALSE,  ignore.case = FALSE,  contrast,  null,  vcov,  coef,  robust = TRUE,  df = NULL,  print = NULL,  labels,  label.width,  only.coef = FALSE,  back.transform = NULL,  folds = 0,  cluster,  R = 0,  null.sim)

Arguments

Details

influence function decomposition of estimator\widehat{\theta} basedon dataZ_1,\ldots,Z_n:

\sqrt{n}(\widehat{\theta}-\theta) =\frac{1}{\sqrt{n}}\sum_{i=1}^n IC(Z_i; P) + o_p(1)

can be extracted withtheIC method.

See Also

estimate.array

Examples

## Simulation from logistic regression modelm <- lvm(y~x+z);distribution(m,y~x) <- binomial.lvm("logit")d <- sim(m,1000)g <- glm(y~z+x,data=d,family=binomial())g0 <- glm(y~1,data=d,family=binomial())## LRTestimate(g,g0)## Plain estimates (robust standard errors)estimate(g)## Testing contrastsestimate(g,null=0)estimate(g,rbind(c(1,1,0),c(1,0,2)))estimate(g,rbind(c(1,1,0),c(1,0,2)),null=c(1,2))estimate(g,2:3) ## same as cbind(0,1,-1)estimate(g,as.list(2:3)) ## same as rbind(c(0,1,0),c(0,0,1))## Alternative syntaxestimate(g,"z","z"-"x",2*"z"-3*"x")estimate(g,"?")  ## Wildcardsestimate(g,"*Int*","z")estimate(g,"1","2"-"3",null=c(0,1))estimate(g,2,3)## Usual (non-robust) confidence intervalsestimate(g,robust=FALSE)## Transformationsestimate(g,function(p) p[1]+p[2])## Multiple parameterse <- estimate(g,function(p) c(p[1]+p[2], p[1]*p[2]))evcov(e)## Label new parametersestimate(g,function(p) list("a1"=p[1]+p[2], "b1"=p[1]*p[2]))##'## Multiple groupm <- lvm(y~x)m <- baptize(m)d2 <- d1 <- sim(m,50,seed=1)e <- estimate(list(m,m),list(d1,d2))estimate(e) ## Wrongee <- estimate(e, id=rep(seq(nrow(d1)), 2)) ## Clusteredeeestimate(lm(y~x,d1))## Marginalizef <- function(p,data)  list(p0=lava:::expit(p["(Intercept)"] + p["z"]*data[,"z"]),       p1=lava:::expit(p["(Intercept)"] + p["x"] + p["z"]*data[,"z"]))e <- estimate(g, f, average=TRUE)eestimate(e,diff)estimate(e,cbind(1,1))## Clusters and subset (conditional marginal effects)d$id <- rep(seq(nrow(d)/4),each=4)estimate(g,function(p,data)         list(p0=lava:::expit(p[1] + p["z"]*data[,"z"])),         subset=d$z>0, id=d$id, average=TRUE)## More examples with clusters:m <- lvm(c(y1,y2,y3)~u+x)d <- sim(m,10)l1 <- glm(y1~x,data=d)l2 <- glm(y2~x,data=d)l3 <- glm(y3~x,data=d)## Some random id-numbersid1 <- c(1,1,4,1,3,1,2,3,4,5)id2 <- c(1,2,3,4,5,6,7,8,1,1)id3 <- seq(10)## Un-stacked and stacked i.i.d. decompositionIC(estimate(l1,id=id1,stack=FALSE))IC(estimate(l1,id=id1))## Combined i.i.d. decompositione1 <- estimate(l1,id=id1)e2 <- estimate(l2,id=id2)e3 <- estimate(l3,id=id3)(a2 <- merge(e1,e2,e3))## If all models were estimated on the same data we could use the## syntax:## Reduce(merge,estimate(list(l1,l2,l3)))## Same:IC(a1 <- merge(l1,l2,l3,id=list(id1,id2,id3)))IC(merge(l1,l2,l3,id=TRUE)) # one-to-one (same clusters)IC(merge(l1,l2,l3,id=FALSE)) # independence## Monte Carlo approach, simple trend test examplem <- categorical(lvm(),~x,K=5)regression(m,additive=TRUE) <- y~xd <- simulate(m,100,seed=1,'y~x'=0.1)l <- lm(y~-1+factor(x),data=d)f <- function(x) coef(lm(x~seq_along(x)))[2]null <- rep(mean(coef(l)),length(coef(l)))##  just need to make sure we simulate under H0: slope=0estimate(l,f,R=1e2,null.sim=null)estimate(l,f)

Estimation of parameters in a Latent Variable Model (lvm)

Description

Estimate parameters. MLE, IV or user-defined estimator.

Usage

## S3 method for class 'lvm'estimate(  x,  data = parent.frame(),  estimator = NULL,  control = list(),  missing = FALSE,  weights,  weightsname,  data2,  id,  fix,  index = !quick,  graph = FALSE,  messages = lava.options()$messages,  quick = FALSE,  method,  param,  cluster,  p,  ...)

Arguments

Details

A list of parameters controlling the estimation and optimization proceduresis parsed via thecontrol argument. By default Maximum Likelihood isused assuming multivariate normal distributed measurement errors. A listwith one or more of the following elements is expected:

start:

Starting value. The order of the parameters can be shown bycallingcoef (withmean=TRUE) on thelvm-object or withplot(..., labels=TRUE). Note that this requires a check that it isactual the model being estimated, asestimate might add additionalrestriction to the model, e.g. through thefix andexo.fixarguments. Thelvm-object of a fitted model can be extracted with theModel-function.

starterfun:

Starter-function with syntaxfunction(lvm, S, mu). Three builtin functions are available:startvalues,startvalues0,startvalues1, ...

estimator:

String defining which estimator to use (Defaults to“gaussian”)

meanstructure

Logical variable indicatingwhether to fit model with meanstructure.

method:

String pointing toalternative optimizer (e.g.optim to use simulated annealing).

control:

Parameters passed to the optimizer (defaultstats::nlminb).

tol:

Tolerance of optimization constraints on lower limit ofvariance parameters.

Value

Alvmfit-object.

Author(s)

Klaus K. Holst

See Also

estimate.default score, information

Examples

dd <- read.table(header=TRUE,text="x1 x2 x3 0.0 -0.5 -2.5-0.5 -2.0  0.0 1.0  1.5  1.0 0.0  0.5  0.0-2.5 -1.5 -1.0")e <- estimate(lvm(c(x1,x2,x3)~u),dd)## Simulation examplem <- lvm(list(y~v1+v2+v3+v4,c(v1,v2,v3,v4)~x))covariance(m) <- v1~v2+v3+v4dd <- sim(m,10000) ## Simulate 10000 observations from modele <- estimate(m, dd) ## Estimate parameterse## Using just sufficient statisticsn <- nrow(dd)e0 <- estimate(m,data=list(S=cov(dd)*(n-1)/n,mu=colMeans(dd),n=n))rm(dd)## Multiple group analysism <- lvm()regression(m) <- c(y1,y2,y3)~uregression(m) <- u~xd1 <- sim(m,100,p=c("u,u"=1,"u~x"=1))d2 <- sim(m,100,p=c("u,u"=2,"u~x"=-1))mm <- baptize(m)regression(mm,u~x) <- NAcovariance(mm,~u) <- NAintercept(mm,~u) <- NAee <- estimate(list(mm,mm),list(d1,d2))## Missing datad0 <- makemissing(d1,cols=1:2)e0 <- estimate(m,d0,missing=TRUE)e0

Add an observed event time outcome to a latent variable model.

Description

For example, if the model 'm' includes latent event time variablesare called 'T1' and 'T2' and 'C' is the end of follow-up (right censored),then one can specify

Usage

eventTime(object, formula, eventName = "status", ...)

Arguments

Details

eventTime(object=m,formula=ObsTime~min(T1=a,T2=b,C=0,"ObsEvent"))

when data are simulated from the modelone gets 2 new columns:

- "ObsTime": the smallest of T1, T2 and C- "ObsEvent": 'a' if T1 is smallest, 'b' if T2 is smallest and '0' if C is smallest

Note that "ObsEvent" and "ObsTime" are names specified by the user.

Author(s)

Thomas A. Gerds, Klaus K. Holst

Examples

# Right censored survival data without covariatesm0 <- lvm()distribution(m0,"eventtime") <- coxWeibull.lvm(scale=1/100,shape=2)distribution(m0,"censtime") <- coxExponential.lvm(rate=1/10)m0 <- eventTime(m0,time~min(eventtime=1,censtime=0),"status")sim(m0,10)# Alternative specification of the right censored survival outcome## eventTime(m,"Status") <- ~min(eventtime=1,censtime=0)# Cox regression:# lava implements two different parametrizations of the same# Weibull regression model. The first specifies# the effects of covariates as proportional hazard ratios# and works as follows:m <- lvm()distribution(m,"eventtime") <- coxWeibull.lvm(scale=1/100,shape=2)distribution(m,"censtime") <- coxWeibull.lvm(scale=1/100,shape=2)m <- eventTime(m,time~min(eventtime=1,censtime=0),"status")distribution(m,"sex") <- binomial.lvm(p=0.4)distribution(m,"sbp") <- normal.lvm(mean=120,sd=20)regression(m,from="sex",to="eventtime") <- 0.4regression(m,from="sbp",to="eventtime") <- -0.01sim(m,6)# The parameters can be recovered using a Cox regression# routine or a Weibull regression model. E.g.,## Not run:     set.seed(18)    d <- sim(m,1000)    library(survival)    coxph(Surv(time,status)~sex+sbp,data=d)    sr <- survreg(Surv(time,status)~sex+sbp,data=d)    library(SurvRegCensCov)    ConvertWeibull(sr)## End(Not run)# The second parametrization is an accelerated failure time# regression model and uses the function weibull.lvm instead# of coxWeibull.lvm to specify the event time distributions.# Here is an example:ma <- lvm()distribution(ma,"eventtime") <- weibull.lvm(scale=3,shape=1/0.7)distribution(ma,"censtime") <- weibull.lvm(scale=2,shape=1/0.7)ma <- eventTime(ma,time~min(eventtime=1,censtime=0),"status")distribution(ma,"sex") <- binomial.lvm(p=0.4)distribution(ma,"sbp") <- normal.lvm(mean=120,sd=20)regression(ma,from="sex",to="eventtime") <- 0.7regression(ma,from="sbp",to="eventtime") <- -0.008set.seed(17)sim(ma,6)# The regression coefficients of the AFT model# can be tranformed into log(hazard ratios):#  coef.coxWeibull = - coef.weibull / shape.weibull## Not run:     set.seed(17)    da <- sim(ma,1000)    library(survival)    fa <- coxph(Surv(time,status)~sex+sbp,data=da)    coef(fa)    c(0.7,-0.008)/0.7## End(Not run)# The following are equivalent parametrizations# which produce exactly the same random numbers:model.aft <- lvm()distribution(model.aft,"eventtime") <- weibull.lvm(intercept=-log(1/100)/2,sigma=1/2)distribution(model.aft,"censtime") <- weibull.lvm(intercept=-log(1/100)/2,sigma=1/2)sim(model.aft,6,seed=17)model.aft <- lvm()distribution(model.aft,"eventtime") <- weibull.lvm(scale=100^(1/2), shape=2)distribution(model.aft,"censtime") <- weibull.lvm(scale=100^(1/2), shape=2)sim(model.aft,6,seed=17)model.cox <- lvm()distribution(model.cox,"eventtime") <- coxWeibull.lvm(scale=1/100,shape=2)distribution(model.cox,"censtime") <- coxWeibull.lvm(scale=1/100,shape=2)sim(model.cox,6,seed=17)# The minimum of multiple latent times one of them still# being a censoring time, yield# right censored competing risks datamc <- lvm()distribution(mc,~X2) <- binomial.lvm()regression(mc) <- T1~f(X1,-.5)+f(X2,0.3)regression(mc) <- T2~f(X2,0.6)distribution(mc,~T1) <- coxWeibull.lvm(scale=1/100)distribution(mc,~T2) <- coxWeibull.lvm(scale=1/100)distribution(mc,~C) <- coxWeibull.lvm(scale=1/100)mc <- eventTime(mc,time~min(T1=1,T2=2,C=0),"event")sim(mc,6)

fplot

Description

Faster plot via RGL

Usage

fplot(  x,  y,  z = NULL,  xlab,  ylab,  ...,  z.col = topo.colors(64),  data = parent.frame(),  add = FALSE,  aspect = c(1, 1),  zoom = 0.8)

Arguments

Examples

if (interactive()) {data(iris)fplot(Sepal.Length ~ Petal.Length+Species, data=iris, size=2, type="s")}

Read Mplus output

Description

Read Mplus output files

Usage

getMplus(infile = "template.out", coef = TRUE, ...)

Arguments

Author(s)

Klaus K. Holst

See Also

getSAS

Read SAS output

Description

Run SAS code like in the following:

Usage

getSAS(infile, entry = "Parameter Estimates", ...)

Arguments

Details

ODS CSVALL BODY="myest.csv";proc nlmixed data=aj qpoints=2 dampstep=0.5;...run;ODS CSVALL Close;

and read results into R with:

getsas("myest.csv","Parameter Estimates")

Author(s)

Klaus K. Holst

See Also

getMplus

Extract model summaries and GOF statistics for model object

Description

Calculates various GOF statistics for model object including globalchi-squared test statistic and AIC. Extract model-specific mean and variancestructure, residuals and various predicitions.

Usage

gof(object, ...)## S3 method for class 'lvmfit'gof(object, chisq=FALSE, level=0.90, rmsea.threshold=0.05,all=FALSE,...)moments(x,...)## S3 method for class 'lvm'moments(x, p, debug=FALSE, conditional=FALSE, data=NULL, latent=FALSE, ...)## S3 method for class 'lvmfit'logLik(object, p=coef(object),                      data=model.frame(object),                      model=object$estimator,                      weights=Weights(object),                      data2=object$data$data2,                          ...)## S3 method for class 'lvmfit'score(x, data=model.frame(x), p=pars(x), model=x$estimator,                   weights=Weights(x), data2=x$data$data2, ...)## S3 method for class 'lvmfit'information(x,p=pars(x),n=x$data$n,data=model.frame(x),                   model=x$estimator,weights=Weights(x), data2=x$data$data2, ...)

Arguments

Value

Ahtest-object.

Author(s)

Klaus K. Holst

Examples

m <- lvm(list(y~v1+v2+v3+v4,c(v1,v2,v3,v4)~x))set.seed(1)dd <- sim(m,1000)e <- estimate(m, dd)gof(e,all=TRUE,rmsea.threshold=0.05,level=0.9)set.seed(1)m <- lvm(list(c(y1,y2,y3)~u,y1~x)); latent(m) <- ~uregression(m,c(y2,y3)~u) <- "b"d <- sim(m,1000)e <- estimate(m,d)rsq(e)##'rr <- rsq(e,TRUE)rrestimate(rr,contrast=rbind(c(1,-1,0),c(1,0,-1),c(0,1,-1)))

Hubble data

Description

Velocity (v) and distance (D) measures of 36 Type Ia super-novae from theHubble Space Telescope

Format

data.frame

Source

Freedman, W. L., et al. 2001, AstroPhysicalJournal, 553, 47.

Hubble data

Description

Hubble data

Format

data.frame

See Also

hubble

Extract i.i.d. decomposition from model object

Description

This function extracts

Usage

iid(x, ...)

Arguments

Organize several image calls (for visualizing categorical data)

Description

Visualize categorical by group variable

Usage

images(  x,  group,  ncol = 2,  byrow = TRUE,  colorbar = 1,  colorbar.space = 0.1,  label.offset = 0.02,  order = TRUE,  colorbar.border = 0,  main,  rowcol = FALSE,  plotfun = NULL,  axis1,  axis2,  mar,  col = list(c("#EFF3FF", "#BDD7E7", "#6BAED6", "#2171B5"), c("#FEE5D9", "#FCAE91",    "#FB6A4A", "#CB181D"), c("#EDF8E9", "#BAE4B3", "#74C476", "#238B45"), c("#FEEDDE",    "#FDBE85", "#FD8D3C", "#D94701")),  ...)

Arguments

Author(s)

Klaus Holst

Examples

X <- matrix(rbinom(400,3,0.5),20)group <- rep(1:4,each=5)images(X,colorbar=0,zlim=c(0,3))images(X,group=group,zlim=c(0,3))## Not run: images(X,group=group,col=list(RColorBrewer::brewer.pal(4,"Purples"),                               RColorBrewer::brewer.pal(4,"Greys"),                               RColorBrewer::brewer.pal(4,"YlGn"),                               RColorBrewer::brewer.pal(4,"PuBuGn")),colorbar=2,zlim=c(0,3))## End(Not run)images(list(X,X,X,X),group=group,zlim=c(0,3))images(list(X,X,X,X),ncol=1,group=group,zlim=c(0,3))images(list(X,X),group,axis2=c(FALSE,FALSE),axis1=c(FALSE,FALSE),      mar=list(c(0,0,0,0),c(0,0,0,0)),yaxs="i",xaxs="i",zlim=c(0,3))

Data

Description

Description

Format

data.frame

Source

Simulated

Fix mean parameters in 'lvm'-object

Description

Define linear constraints on intercept parameters in alvm-object.

Usage

## S3 replacement method for class 'lvm'intercept(object, vars, ...) <- value

Arguments

Details

Theintercept function is used to specify linear constraints on theintercept parameters of a latent variable model. As an example we look atthe multivariate regression model

E(Y_1|X) = \alpha_1 + \beta_1 X

E(Y_2|X) = \alpha_2 + \beta_2X

defined by the call

m <- lvm(c(y1,y2) ~ x)

To fix\alpha_1=\alpha_2 we call

intercept(m) <- c(y1,y2) ~ f(mu)

Fixed parameters can be reset by fixing them toNA. For instance tofree the parameter restriction ofY_1 and at the same time fixing\alpha_2=2, we call

intercept(m, ~y1+y2) <- list(NA,2)

Callingintercept with no additional arguments will return thecurrent intercept restrictions of thelvm-object.

Value

Alvm-object

Note

Variables will be added to the model if not already present.

Author(s)

Klaus K. Holst

See Also

covariance<-,regression<-,constrain<-,parameter<-,latent<-,cancel<-,kill<-

Examples

## A multivariate modelm <- lvm(c(y1,y2) ~ f(x1,beta)+x2)regression(m) <- y3 ~ f(x1,beta)intercept(m) <- y1 ~ f(mu)intercept(m, ~y2+y3) <- list(2,"mu")intercept(m) ## Examine intercepts of model (NA translates to free/unique paramete##r)

Define intervention

Description

Define intervention in a 'lvm' object

Usage

## S3 method for class 'lvm'intervention(object, to, value, dist = none.lvm(), ...)

Arguments

See Also

regression lvm sim

Examples

m <- lvm(y ~ a + x, a ~ x)distribution(m, ~a+y) <- binomial.lvm()mm <- intervention(m, "a", value=3)sim(mm, 10)mm <- intervention(m, a~x, function(x) (x>0)*1)sim(mm, 10)

Plot/estimate surface

Description

Plot/estimate surface

Usage

ksmooth2(  x,  data,  h = NULL,  xlab = NULL,  ylab = NULL,  zlab = "",  gridsize = rep(51L, 2),  ...)

Arguments

Examples

if (requireNamespace("KernSmooth")) {##'ksmooth2(rmvn0(1e4,sigma=diag(2)*.5+.5),c(-3.5,3.5),h=1,        rgl=FALSE,theta=30)##'if (interactive()) {    ksmooth2(rmvn0(1e4,sigma=diag(2)*.5+.5),c(-3.5,3.5),h=1)    ksmooth2(function(x,y) x^2+y^2, c(-20,20))    ksmooth2(function(x,y) x^2+y^2, xlim=c(-5,5), ylim=c(0,10))    f <- function(x,y) 1-sqrt(x^2+y^2)    surface(f,xlim=c(-1,1),alpha=0.9,aspect=c(1,1,0.75))    surface(f,xlim=c(-1,1),clut=heat.colors(128))    ##play3d(spin3d(axis=c(0,0,1), rpm=8), duration=5)}if (interactive()) {    surface(function(x) dmvn0(x,sigma=diag(2)),c(-3,3),lit=FALSE,smooth=FALSE,box=FALSE,alpha=0.8)    surface(function(x) dmvn0(x,sigma=diag(2)),c(-3,3),box=FALSE,specular="black")##'}if (!inherits(try(find.package("fields"),silent=TRUE),"try-error")) {    f <- function(x,y) 1-sqrt(x^2+y^2)    ksmooth2(f,c(-1,1),rgl=FALSE,image=fields::image.plot)}}

Define labels of graph

Description

Alters labels of nodes and edges in the graph of a latent variable model

Usage

## Default S3 replacement method:labels(object, ...) <- value## S3 replacement method for class 'lvm'edgelabels(object, to, ...) <- value## Default S3 replacement method:nodecolor(object, var=vars(object),border, labcol, shape, lwd, ...) <- value

Arguments

Author(s)

Klaus K. Holst

Examples

m <- lvm(c(y,v)~x+z)regression(m) <- c(v,x)~zlabels(m) <- c(y=expression(psi), z=expression(zeta))nodecolor(m,~y+z+x,border=c("white","white","black"),          labcol="white", lwd=c(1,1,5),          lty=c(1,2)) <-  c("orange","indianred","lightgreen")edgelabels(m,y~z+x, cex=c(2,1.5), col=c("orange","black"),labcol="darkblue",           arrowhead=c("tee","dot"),           lwd=c(3,1)) <- expression(phi,rho)edgelabels(m,c(v,x)~z, labcol="red", cex=0.8,arrowhead="none") <- 2if (interactive()) {    plot(m,addstyle=FALSE)}m <- lvm(y~x)labels(m) <- list(x="multiple\nlines")if (interactive()) {op <- par(mfrow=c(1,2))plot(m,plain=TRUE)plot(m)par(op)d <- sim(m,100)e <- estimate(m,d)plot(e,type="sd")}

Set global options forlava

Description

Extract and set global parameters oflava. In particular optimizationparameters for theestimate function.

Usage

lava.options(...)

Arguments

Details

• param: 'relative' (factor loading and variance of oneendogenous variables in each measurement model are fixed to one), 'absolute'(mean and variance of latent variables are set to 0 and 1, respectively),'hybrid' (intercept of latent variables is fixed to 0, and factor loading ofat least one endogenous variable in each measurement model is fixed to 1),'none' (no constraints are added)

• layout: One of 'dot','fdp','circo','twopi','neato','osage'

• messages: Set to 0 to disable various output messages

• ...

seecontrol parameter of theestimate function.

Value

list of parameters

Author(s)

Klaus K. Holst

Examples

## Not run: lava.options(iter.max=100,messages=0)## End(Not run)

Initialize new latent variable model

Description

Function that constructs a new latent variable model object

Usage

lvm(x = NULL, ..., latent = NULL, messages = lava.options()$messages)

Arguments

Value

Returns an object of classlvm.

Author(s)

Klaus K. Holst

See Also

regression,covariance,intercept, ...

Examples

m <- lvm() # Empty modelm1 <- lvm(y~x) # Simple linear regressionm2 <- lvm(~y1+y2) # Model with two independent variables (argument)m3 <- lvm(list(c(y1,y2,y3)~u,u~x+z)) # SEM with three items

Create random missing data

Description

Generates missing entries in data.frame/matrix

Usage

makemissing(  data,  p = 0.2,  cols = seq_len(ncol(data)),  rowwise = FALSE,  nafun = function(x) x,  seed = NULL)

Arguments

Value

data.frame

Author(s)

Klaus K. Holst

Two-stage (non-linear) measurement error

Description

Two-stage measurement error

Usage

measurement.error(  model1,  formula,  data = parent.frame(),  predictfun = function(mu, var, data, ...) mu[, 1]^2 + var[1],  id1,  id2,  ...)

Arguments

See Also

stack.estimate

Examples

m <- lvm(c(y1,y2,y3)~u,c(y3,y4,y5)~v,u~~v,c(u,v)~x)transform(m,u2~u) <- function(x) x^2transform(m,uv~u+v) <- prodregression(m) <- z~u2+u+v+uv+xset.seed(1)d <- sim(m,1000,p=c("u,u"=1))## Stage 1m1 <- lvm(c(y1[0:s],y2[0:s],y3[0:s])~1*u,c(y3[0:s],y4[0:s],y5[0:s])~1*v,u~b*x,u~~v)latent(m1) <- ~u+ve1 <- estimate(m1,d)pp <- function(mu,var,data,...) {    cbind(u=mu[,"u"],u2=mu[,"u"]^2+var["u","u"],v=mu[,"v"],uv=mu[,"u"]*mu[,"v"]+var["u","v"])}(e <- measurement.error(e1, z~1+x, data=d, predictfun=pp))## uu <- seq(-1,1,length.out=100)## pp <- estimate(e,function(p,...) p["(Intercept)"]+p["u"]*uu+p["u2"]*uu^2)$coefmatif (interactive()) {    plot(e,intercept=TRUE,line=0)    f <- function(p) p[1]+p["u"]*u+p["u2"]*u^2    u <- seq(-1,1,length.out=100)    plot(e, f, data=data.frame(u), ylim=c(-.5,2.5))}

Missing data example

Description

Simulated data generated from model

E(Y_i\mid X) = X, \quad cov(Y_1,Y_2\mid X)=0.5

Format

list of data.frames

Details

The list contains four data sets1) Complete data2) MCAR3) MAR4) MNAR (missing mechanism depends on variable V correlated with Y1,Y2)

Source

Simulated

Examples

data(missingdata)e0 <- estimate(lvm(c(y1,y2)~b*x,y1~~y2),missingdata[[1]]) ## No missinge1 <- estimate(lvm(c(y1,y2)~b*x,y1~~y2),missingdata[[2]]) ## CC (MCAR)e2 <- estimate(lvm(c(y1,y2)~b*x,y1~~y2),missingdata[[2]],missing=TRUE) ## MCARe3 <- estimate(lvm(c(y1,y2)~b*x,y1~~y2),missingdata[[3]]) ## CC (MAR)e4 <- estimate(lvm(c(y1,y2)~b*x,y1~~y2),missingdata[[3]],missing=TRUE) ## MAR

Estimate mixture latent variable model.

Description

Estimate mixture latent variable model

Usage

mixture(  x,  data,  k = length(x),  control = list(),  vcov = "observed",  names = FALSE,  ...)

Arguments

Details

Estimate parameters in a mixture of latent variable models via the EMalgorithm.

The performance of the EM algorithm can be tuned via thecontrolargument, a list where a subset of the following members can be altered:

start

Optional starting values

nstart

Evaluatenstart different starting values and run the EM-algorithm on theparameters with largest likelihood

tol

Convergence tolerance of theEM-algorithm. The algorithm is stopped when the absolute change inlikelihood and parameter (2-norm) between successive iterations is less thantol

iter.max

Maximum number of iterations of theEM-algorithm

gamma

Scale-down (i.e. number between 0 and 1) of thestep-size of the Newton-Raphson algorithm in the M-step

trace

Traceinformation on the EM-algorithm is printed on everytracethiteration

Note that the algorithm can be aborted any time (C-c) and still be saved(via on.exit call).

Author(s)

Klaus K. Holst

See Also

mvnmix

Examples

m0 <- lvm(list(y~x+z,x~z))distribution(m0,~z) <- binomial.lvm()d <- sim(m0,2000,p=c("y~z"=2,"y~x"=1),seed=1)## unmeasured confounder examplem <- baptize(lvm(y~x, x~1));intercept(m,~x+y) <- NAif (requireNamespace('mets', quietly=TRUE)) {  set.seed(42)  M <- mixture(m,k=2,data=d,control=list(trace=1,tol=1e-6))  summary(M)  lm(y~x,d)  estimate(M,"y~x")  ## True slope := 1}

Model searching

Description

Performs Wald or score tests

Usage

modelsearch(x, k = 1, dir = "forward", type = "all", ...)

Arguments

Value

Matrix of test-statistics and p-values

Author(s)

Klaus K. Holst

See Also

compare,equivalence

Examples

m <- lvm();regression(m) <- c(y1,y2,y3) ~ eta; latent(m) <- ~etaregression(m) <- eta ~ xm0 <- m; regression(m0) <- y2 ~ xdd <- sim(m0,100)[,manifest(m0)]e <- estimate(m,dd);modelsearch(e,messages=0)modelsearch(e,messages=0,type="cor")

Estimate probabilities in contingency table

Description

Estimate probabilities in contingency table

Usage

multinomial(  x,  data = parent.frame(),  marginal = FALSE,  transform,  vcov = TRUE,  IC = TRUE,  ...)

Arguments

Author(s)

Klaus K. Holst

Examples

set.seed(1)breaks <- c(-Inf,-1,0,Inf)m <- lvm(); covariance(m,pairwise=TRUE) <- ~y1+y2+y3+y4d <- transform(sim(m,5e2),              z1=cut(y1,breaks=breaks),              z2=cut(y2,breaks=breaks),              z3=cut(y3,breaks=breaks),              z4=cut(y4,breaks=breaks))multinomial(d[,5])(a1 <- multinomial(d[,5:6]))(K1 <- kappa(a1)) ## Cohen's kappaK2 <- kappa(d[,7:8])## Testing difference K1-K2:estimate(merge(K1,K2,id=TRUE),diff)estimate(merge(K1,K2,id=FALSE),diff) ## Wrong std.err ignoring dependencesqrt(vcov(K1)+vcov(K2))## Average of the two kappas:estimate(merge(K1,K2,id=TRUE),function(x) mean(x))estimate(merge(K1,K2,id=FALSE),function(x) mean(x)) ## Independence##'## Goodman-Kruskal's gammam2 <- lvm(); covariance(m2) <- y1~y2breaks1 <- c(-Inf,-1,0,Inf)breaks2 <- c(-Inf,0,Inf)d2 <- transform(sim(m2,5e2),              z1=cut(y1,breaks=breaks1),              z2=cut(y2,breaks=breaks2))(g1 <- gkgamma(d2[,3:4]))## same as## Not run: gkgamma(table(d2[,3:4]))gkgamma(multinomial(d2[,3:4]))## End(Not run)##partial gammad2$x <- rbinom(nrow(d2),2,0.5)gkgamma(z1~z2|x,data=d2)

Estimate mixture latent variable model

Description

Estimate mixture latent variable model

Usage

mvnmix(  data,  k = 2,  theta,  steps = 500,  tol = 1e-16,  lambda = 0,  mu = NULL,  silent = TRUE,  extra = FALSE,  n.start = 1,  init = "kmpp",  ...)

Arguments

Details

Estimate parameters in a mixture of latent variable models via the EMalgorithm.

Value

Amixture object

Author(s)

Klaus K. Holst

See Also

mixture

Examples

data(faithful)set.seed(1)M1 <- mvnmix(faithful[,"waiting",drop=FALSE],k=2)M2 <- mvnmix(faithful,k=2)if (interactive()) {    par(mfrow=c(2,1))    plot(M1,col=c("orange","blue"),ylim=c(0,0.05))    plot(M2,col=c("orange","blue"))}

Example data (nonlinear model)

Description

Example data (nonlinear model)

Format

data.frame

Source

Simulated

Example SEM data (nonlinear)

Description

Simulated data

Format

data.frame

Source

Simulated

Define variables as ordinal

Description

Define variables as ordinal in latent variable model object

Usage

ordinal(x, ...) <- value

Arguments

Examples

if (requireNamespace("mets")) {m <- lvm(y + z ~ x + 1*u[0], latent=~u)ordinal(m, K=3) <- ~y+zd <- sim(m, 100, seed=1)e <- estimate(m, d)}

Univariate cumulative link regression models

Description

Ordinal regression models

Usage

ordreg(  formula,  data = parent.frame(),  offset,  family = stats::binomial("probit"),  start,  fast = FALSE,  ...)

Arguments

Author(s)

Klaus K. Holst

Examples

m <- lvm(y~x)ordinal(m,K=3) <- ~yd <- sim(m,100)e <- ordreg(y~x,d)

Generic method for finding indeces of model parameters

Description

Generic method for finding indeces of model parameters

Usage

parpos(x, ...)

Arguments

Author(s)

Klaus K. Holst

Calculate partial correlations

Description

Calculate partial correlation coefficients and confidence limits via Fishersz-transform

Usage

partialcor(formula, data, level = 0.95, ...)

Arguments

Value

A coefficient matrix

Author(s)

Klaus K. Holst

Examples

m <- lvm(c(y1,y2,y3)~x1+x2)covariance(m) <- c(y1,y2,y3)~y1+y2+y3d <- sim(m,500)partialcor(~x1+x2,d)

Extract pathways in model graph

Description

Extract all possible paths from one variable to another connected componentin a latent variable model. In an estimated model the effect size isdecomposed into direct, indirect and total effects including approximatestandard errors.

Usage

## S3 method for class 'lvm' path(object, to = NULL, from, all=FALSE, ...)## S3 method for class 'lvmfit' effects(object, to, from, ...)

Arguments

Value

Ifobject is of classlvmfit a list with the followingelements is returned

Ifobject is of classlvm only thepath element will bereturned.

Theeffects method returns an object of classeffects.

Note

For alvmfit-object the parameters estimates and theircorresponding covariance matrix are also returned. Theeffects-function additionally calculates the total and indirecteffects with approximate standard errors

Author(s)

Klaus K. Holst

See Also

children,parents

Examples

m <- lvm(c(y1,y2,y3)~eta)regression(m) <- y2~x1latent(m) <- ~etaregression(m) <- eta~x1+x2d <- sim(m,500)e <- estimate(m,d)path(Model(e),y2~x1)parents(Model(e), ~y2)children(Model(e), ~x2)children(Model(e), ~x2+eta)effects(e,y2~x1)## All simple paths (undirected)path(m,y1~x1,all=TRUE)

Polychoric correlation

Description

Maximum likelhood estimates of polychoric correlations

Usage

pcor(x, y, X, start, ...)

Arguments

Convert pdf to raster format

Description

Convert PDF file to print quality png (default 300 dpi)

Usage

pdfconvert(  files,  dpi = 300,  resolution = 1024,  gs,  gsopt,  resize,  format = "png",  ...)

Arguments

Details

Access to ghostscript program 'gs' is needed

Author(s)

Klaus K. Holst

See Also

dev.copy2pdf,printdev

Plot method for 'estimate' objects

Description

Plot method for 'estimate' objects

Usage

## S3 method for class 'estimate'plot(  x,  f,  idx,  intercept = FALSE,  data,  confint = TRUE,  type = "l",  xlab = "x",  ylab = "f(x)",  col = 1,  add = FALSE,  ...)

Arguments

Plot path diagram

Description

Plot the path diagram of a SEM

Usage

## S3 method for class 'lvm'plot(  x,  diag = FALSE,  cor = TRUE,  labels = FALSE,  intercept = FALSE,  addcolor = TRUE,  plain = FALSE,  cex,  fontsize1 = 10,  noplot = FALSE,  graph = list(rankdir = "BT"),  attrs = list(graph = graph),  unexpr = FALSE,  addstyle = TRUE,  plot.engine = lava.options()$plot.engine,  init = TRUE,  layout = lava.options()$layout,  edgecolor = lava.options()$edgecolor,  graph.proc = lava.options()$graph.proc,  ...)

Arguments

Author(s)

Klaus K. Holst

Examples

if (interactive()) {m <- lvm(c(y1,y2) ~ eta)regression(m) <- eta ~ z+x2regression(m) <- c(eta,z) ~ x1latent(m) <- ~etalabels(m) <- c(y1=expression(y[scriptscriptstyle(1)]),y2=expression(y[scriptscriptstyle(2)]),x1=expression(x[scriptscriptstyle(1)]),x2=expression(x[scriptscriptstyle(2)]),eta=expression(eta))edgelabels(m, eta ~ z+x1+x2, cex=2, lwd=3,           col=c("orange","lightblue","lightblue")) <- expression(rho,phi,psi)nodecolor(m, vars(m), border="white", labcol="darkblue") <- NAnodecolor(m, ~y1+y2+z, labcol=c("white","white","black")) <- NAplot(m,cex=1.5)d <- sim(m,100)e <- estimate(m,d)plot(e)m <- lvm(c(y1,y2) ~ eta)regression(m) <- eta ~ z+x2regression(m) <- c(eta,z) ~ x1latent(m) <- ~etaplot(lava:::beautify(m,edgecol=FALSE))}

Plot method for simulation 'sim' objects

Description

Density and scatter plots

Usage

## S3 method for class 'sim'plot(  x,  estimate,  se = NULL,  true = NULL,  names = NULL,  auto.layout = TRUE,  byrow = FALSE,  type = "p",  ask = grDevices::dev.interactive(),  col = c("gray60", "orange", "darkblue", "seagreen", "darkred"),  pch = 16,  cex = 0.5,  lty = 1,  lwd = 0.3,  legend,  legendpos = "topleft",  cex.legend = 0.8,  plot.type = c("multiple", "single"),  polygon = TRUE,  density = 0,  angle = -45,  cex.axis = 0.8,  alpha = 0.2,  main,  cex.main = 1,  equal = FALSE,  delta = 1.15,  ylim = NULL,  xlim = NULL,  ylab = "",  xlab = "",  rug = FALSE,  rug.alpha = 0.5,  line.col = scatter.col,  line.lwd = 1,  line.lty = 1,  line.alpha = 1,  scatter.ylab = "Estimate",  scatter.ylim = NULL,  scatter.xlim = NULL,  scatter.alpha = 0.5,  scatter.col = col,  border = col,  true.lty = 2,  true.col = "gray70",  true.lwd = 1.2,  density.plot = TRUE,  scatter.plot = FALSE,  running.mean = scatter.plot,  ...)

Arguments

Examples

n <- 1000val <- cbind(est1=rnorm(n,sd=1),est2=rnorm(n,sd=0.2),est3=rnorm(n,1,sd=0.5),             sd1=runif(n,0.8,1.2),sd2=runif(n,0.1,0.3),sd3=runif(n,0.25,0.75))plot.sim(val,estimate=c(1,2),true=c(0,0),se=c(4,5),equal=TRUE,scatter.plot=TRUE)plot.sim(val,estimate=c(1,3),true=c(0,1),se=c(4,6),xlim=c(-3,3),scatter.ylim=c(-3,3),scatter.plot=TRUE)plot.sim(val,estimate=c(1,2),true=c(0,0),se=c(4,5),equal=TRUE,plot.type="single",scatter.plot=TRUE)plot.sim(val,estimate=c(1),se=c(4,5,6),plot.type="single",scatter.plot=TRUE)plot.sim(val,estimate=c(1,2,3),equal=TRUE,scatter.plot=TRUE)plot.sim(val,estimate=c(1,2,3),equal=TRUE,byrow=TRUE,scatter.plot=TRUE)plot.sim(val,estimate=c(1,2,3),plot.type="single",scatter.plot=TRUE)plot.sim(val,estimate=1,se=c(3,4,5),plot.type="single",scatter.plot=TRUE)density.sim(val,estimate=c(1,2,3),density=c(0,10,10), lwd=2, angle=c(0,45,-45),cex.legend=1.3)

Plot regression lines

Description

Plot regression line (with interactions) and partial residuals.

Usage

plotConf(  model,  var1 = NULL,  var2 = NULL,  data = NULL,  ci.lty = 0,  ci = TRUE,  level = 0.95,  pch = 16,  lty = 1,  lwd = 2,  npoints = 100,  xlim,  col = NULL,  colpt,  alpha = 0.5,  cex = 1,  delta = 0.07,  centermark = 0.03,  jitter = 0.2,  cidiff = FALSE,  mean = TRUE,  legend = ifelse(is.null(var1), FALSE, "topright"),  trans = function(x) {     x },  partres = inherits(model, "lm"),  partse = FALSE,  labels,  vcov,  predictfun,  plot = TRUE,  new = TRUE,  ...)

Arguments

Value

list with following members:

Author(s)

Klaus K. Holst

See Also

termplot

Examples

n <- 100x0 <- rnorm(n)x1 <- seq(-3,3, length.out=n)x2 <- factor(rep(c(1,2),each=n/2), labels=c("A","B"))y <- 5 + 2*x0 + 0.5*x1 + -1*(x2=="B")*x1 + 0.5*(x2=="B") + rnorm(n, sd=0.25)dd <- data.frame(y=y, x1=x1, x2=x2)lm0 <- lm(y ~ x0 + x1*x2, dd)plotConf(lm0, var1="x1", var2="x2")abline(a=5,b=0.5,col="red")abline(a=5.5,b=-0.5,col="red")### points(5+0.5*x1 -1*(x2=="B")*x1 + 0.5*(x2=="B") ~ x1, cex=2)data(iris)l <- lm(Sepal.Length ~ Sepal.Width*Species,iris)plotConf(l,var2="Species")plotConf(l,var1="Sepal.Width",var2="Species")## Not run: ## lme4 modeldd$Id <- rbinom(n, size = 3, prob = 0.3)lmer0 <- lme4::lmer(y ~ x0 + x1*x2 + (1|Id), dd)plotConf(lmer0, var1="x1", var2="x2")## End(Not run)

Prediction in structural equation models

Description

Prediction in structural equation models

Usage

## S3 method for class 'lvm'predict(  object,  x = NULL,  y = NULL,  residual = FALSE,  p,  data,  path = FALSE,  quick = is.null(x) & !(residual | path),  ...)

Arguments

See Also

predictlvm

Examples

m <- lvm(list(c(y1,y2,y3)~u,u~x)); latent(m) <- ~ud <- sim(m,100)e <- estimate(m,d)## Conditional mean (and variance as attribute) given covariatesr <- predict(e)## Best linear unbiased predictor (BLUP)r <- predict(e,vars(e))##  Conditional mean of y3 giving covariates and y1,y2r <- predict(e,y3~y1+y2)##  Conditional mean  gives covariates and y1r <- predict(e,~y1)##  Predicted residuals (conditional on all observed variables)r <- predict(e,vars(e),residual=TRUE)

Predict function for latent variable models

Description

Predictions of conditinoal mean and variance and calculation ofjacobian with respect to parameter vector.

Usage

predictlvm(object, formula, p = coef(object), data = model.frame(object), ...)

Arguments

See Also

predict.lvm

Examples

m <- lvm(c(x1,x2,x3)~u1,u1~z,         c(y1,y2,y3)~u2,u2~u1+z)latent(m) <- ~u1+u2d <- simulate(m,10,"u2,u2"=2,"u1,u1"=0.5,seed=123)e <- estimate(m,d)## Conditional mean given covariatespredictlvm(e,c(x1,x2)~1)$mean## Conditional variance of u1,y1 given x1,x2predictlvm(e,c(u1,y1)~x1+x2)$var

AppendingSurv objects

Description

rbind method forSurv objects

Usage

## S3 method for class 'Surv'rbind(...)

Arguments

Value

Surv object

Author(s)

Klaus K. Holst

Examples

y <- yl <- yr <- rnorm(10)yl[1:5] <- NA; yr[6:10] <- NAS1 <- survival::Surv(yl,yr,type="interval2")S2 <- survival::Surv(y,y>0,type="right")S3 <- survival::Surv(y,y<0,type="left")rbind(S1,S1)rbind(S2,S2)rbind(S3,S3)

Add regression association to latent variable model

Description

Define regression association between variables in alvm-object anddefine linear constraints between model equations.

Usage

## S3 method for class 'lvm'regression(object = lvm(), to, from, fn = NA,messages = lava.options()$messages, additive=TRUE, y, x, value, ...)## S3 replacement method for class 'lvm'regression(object, to=NULL, quick=FALSE, ...) <- value

Arguments

Details

Theregression function is used to specify linear associationsbetween variables of a latent variable model, and offers formula syntaxresembling the model specification of e.g.lm.

For instance, to add the following linear regression model, to thelvm-object,m:

E(Y|X_1,X_2) = \beta_1 X_1 + \beta_2 X_2

We can write

regression(m) <- y ~ x1 + x2

Multivariate models can be specified by successive calls withregression, but multivariate formulas are also supported, e.g.

regression(m) <- c(y1,y2) ~ x1 + x2

defines

E(Y_i|X_1,X_2) = \beta_{1i} X_1 + \beta_{2i} X_2

The special function,f, can be used in the model specification tospecify linear constraints. E.g. to fix\beta_1=\beta_2, we could write

regression(m) <- y ~ f(x1,beta) + f(x2,beta)

The second argument off can also be a number (e.g. defining anoffset) or be set toNA in order to clear any previously definedlinear constraints.

Alternatively, a more straight forward notation can be used:

regression(m) <- y ~ beta*x1 + beta*x2

All the parameter values of the linear constraints can be given as the righthandside expression of the assigment functionregression<- (orregfix<-) if the first (and possibly second) argument is defined aswell. E.g:

regression(m,y1~x1+x2) <- list("a1","b1")

definesE(Y_1|X_1,X_2) = a1 X_1 + b1 X_2. The rhs argument can be amixture of character and numeric values (and NA's to remove constraints).

The functionregression (called without additional arguments) can beused to inspect the linear constraints of alvm-object.

Value

Alvm-object

Note

Variables will be added to the model if not already present.

Author(s)

Klaus K. Holst

See Also

intercept<-,covariance<-,constrain<-,parameter<-,latent<-,cancel<-,kill<-

Examples

m <- lvm() ## Initialize empty lvm-object### E(y1|z,v) = beta1*z + beta2*vregression(m) <- y1 ~ z + v### E(y2|x,z,v) = beta*x + beta*z + 2*v + beta3*uregression(m) <- y2 ~ f(x,beta) + f(z,beta)  + f(v,2) + u### Clear restriction on association between y and### fix slope coefficient of u to betaregression(m, y2 ~ v+u) <- list(NA,"beta")regression(m) ## Examine current linear parameter constraints## ## A multivariate model, E(yi|x1,x2) = beta[1i]*x1 + beta[2i]*x2:m2 <- lvm(c(y1,y2) ~ x1+x2)

Create/extract 'reverse'-diagonal matrix or off-diagonal elements

Description

Create/extract 'reverse'-diagonal matrix or off-diagonal elements

Usage

revdiag(x,...)offdiag(x,type=0,...)revdiag(x,...) <- valueoffdiag(x,type=0,...) <- value

Arguments

Author(s)

Klaus K. Holst

Remove variables from (model) object.

Description

Generic method for removing elements of object

Usage

rmvar(x, ...) <- value

Arguments

Author(s)

Klaus K. Holst

See Also

cancel

Examples

m <- lvm()addvar(m) <- ~y1+y2+xcovariance(m) <- y1~y2regression(m) <- c(y1,y2) ~ x### Cancel the covariance between the residuals of y1 and y2cancel(m) <- y1~y2### Remove y2 from the modelrmvar(m) <- ~y2

Performs a rotation in the plane

Description

Performs a rotation in the plane

Usage

rotate2(x, theta = pi)

Arguments

Value

Returns a matrix of the same dimension asx

Author(s)

Klaus K. Holst

Examples

rotate2(cbind(c(1,2),c(2,1)))

Calculate simultaneous confidence limits by Scheffe's method

Description

Function to compute the Scheffe corrected confidenceinterval for the regression line

Usage

scheffe(model, newdata = model.frame(model), level = 0.95)

Arguments

Examples

x <- rnorm(100)d <- data.frame(y=rnorm(length(x),x),x=x)l <- lm(y~x,d)plot(y~x,d)abline(l)d0 <- data.frame(x=seq(-5,5,length.out=100))d1 <- cbind(d0,predict(l,newdata=d0,interval="confidence"))d2 <- cbind(d0,scheffe(l,d0))lines(lwr~x,d1,lty=2,col="red")lines(upr~x,d1,lty=2,col="red")lines(lwr~x,d2,lty=2,col="blue")lines(upr~x,d2,lty=2,col="blue")

Example SEM data

Description

Simulated data

Format

data.frame

Source

Simulated

Serotonin data

Description

This simulated data mimics a PET imaging study where the 5-HT2Areceptor and serotonin transporter (SERT) binding potential hasbeen quantified into 8 different regions. The 5-HT2Acortical regions are considered high-binding regionsmeasurements. These measurements can be regarded as proxy measures ofthe extra-cellular levels of serotonin in the brain

Format

data.frame

Source

Simulated

Monte Carlo simulation

Description

Applies a function repeatedly for a specified number of replications or overa list/data.frame with plot and summary methods for summarizing the MonteCarlo experiment. Can be parallelized via the future package (use thefuture::plan function).

Usage

## Default S3 method:sim(  x = NULL,  R = 100,  f = NULL,  colnames = NULL,  seed = NULL,  args = list(),  iter = FALSE,  mc.cores,  progressr.message = NULL,  ...)

Arguments

Details

To parallelize the calculation use the future::plan function (e.g.,future::plan(multisession()) to distribute the calculations over the Rreplications on all available cores). The output is controlled via theprogressr package (e.g., progressr::handlers(global=TRUE) to enableprogress information).

See Also

summary.sim plot.sim print.sim sim.lvm

Examples

m <- lvm(y~x+e)distribution(m,~y) <- 0distribution(m,~x) <- uniform.lvm(a=-1.1,b=1.1)transform(m,e~x) <- function(x) (1*x^4)*rnorm(length(x),sd=1)onerun <- function(iter=NULL,...,n=2e3,b0=1,idx=2) {    d <- sim(m,n,p=c("y~x"=b0))    l <- lm(y~x,d)    res <- c(coef(summary(l))[idx,1:2],             confint(l)[idx,],             estimate(l,only.coef=TRUE)[idx,2:4])    names(res) <- c("Estimate","Model.se","Model.lo","Model.hi",                    "Sandwich.se","Sandwich.lo","Sandwich.hi")    res}val <- sim(onerun,R=10,b0=1)valval <- sim(val,R=40,b0=1) ## append resultssummary(val,estimate=c(1,1),confint=c(3,4,6,7),true=c(1,1))summary(val,estimate=c(1,1),se=c(2,5),names=c("Model","Sandwich"))summary(val,estimate=c(1,1),se=c(2,5),true=c(1,1),        names=c("Model","Sandwich"),confint=TRUE)if (interactive()) {    plot(val,estimate=1,c(2,5),true=1,         names=c("Model","Sandwich"),polygon=FALSE)    plot(val,estimate=c(1,1),se=c(2,5),main=NULL,         true=c(1,1),names=c("Model","Sandwich"),         line.lwd=1,col=c("gray20","gray60"),         rug=FALSE)    plot(val,estimate=c(1,1),se=c(2,5),true=c(1,1),         names=c("Model","Sandwich"))}f <- function(a=1, b=1) {  rep(a*b, 5)}R <- Expand(a=1:3, b=1:3)sim(f, R)sim(function(a,b) f(a,b), 3, args=c(a=5,b=5))sim(function(iter=1,a=5,b=5) iter*f(a,b), iter=TRUE, R=5)

Simulate model

Description

Simulate data from a general SEM model including non-linear effects andgeneral link and distribution of variables.

Usage

## S3 method for class 'lvm'sim(x, n = NULL, p = NULL, normal = FALSE, cond = FALSE,sigma = 1, rho = 0.5, X = NULL, unlink=FALSE, latent=TRUE,use.labels = TRUE, seed=NULL, ...)

Arguments

Author(s)

Klaus K. Holst

Examples

#################################################### Logistic regression##################################################m <- lvm(y~x+z)regression(m) <- x~zdistribution(m,~y+z) <- binomial.lvm("logit")d <- sim(m,1e3)head(d)e <- estimate(m,d,estimator="glm")e## Simulate a few observation from estimated modelsim(e,n=5)#################################################### Poisson##################################################distribution(m,~y) <- poisson.lvm()d <- sim(m,1e4,p=c(y=-1,"y~x"=2,z=1))head(d)estimate(m,d,estimator="glm")mean(d$z); lava:::expit(1)summary(lm(y~x,sim(lvm(y[1:2]~4*x),1e3)))##################################################### Gamma distribution##################################################m <- lvm(y~x)distribution(m,~y+x) <- list(Gamma.lvm(shape=2),binomial.lvm())intercept(m,~y) <- 0.5d <- sim(m,1e4)summary(g <- glm(y~x,family=Gamma(),data=d))## Not run: MASS::gamma.shape(g)args(lava::Gamma.lvm)distribution(m,~y) <- Gamma.lvm(shape=2,log=TRUE)sim(m,10,p=c(y=0.5))[,"y"]##################################################### Beta##################################################m <- lvm()distribution(m,~y) <- beta.lvm(alpha=2,beta=1)var(sim(m,100,"y,y"=2))distribution(m,~y) <- beta.lvm(alpha=2,beta=1,scale=FALSE)var(sim(m,100))##################################################### Transform##################################################m <- lvm()transform(m,xz~x+z) <- function(x) x[1]*(x[2]>0)regression(m) <- y~x+z+xzd <- sim(m,1e3)summary(lm(y~x+z + x*I(z>0),d))##################################################### Non-random variables##################################################m <- lvm()distribution(m,~x+z+v+w) <- list(Sequence.lvm(0,5),## Seq. 0 to 5 by 1/n                               Binary.lvm(),       ## Vector of ones                               Binary.lvm(0.5),    ##  0.5n 0, 0.5n 1                               Binary.lvm(interval=list(c(0.3,0.5),c(0.8,1))))sim(m,10)##################################################### Cox model### piecewise constant hazard################################################m <- lvm(t~x)rates <- c(1,0.5); cuts <- c(0,5)## Constant rate: 1 in [0,5), 0.5 in [5,Inf)distribution(m,~t) <- coxExponential.lvm(rate=rates,timecut=cuts)## Not run:     d <- sim(m,2e4,p=c("t~x"=0.1)); d$status <- TRUE    plot(timereg::aalen(survival::Surv(t,status)~x,data=d,                        resample.iid=0,robust=0),spec=1)    L <- approxfun(c(cuts,max(d$t)),f=1,                   cumsum(c(0,rates*diff(c(cuts,max(d$t))))),                   method="linear")    curve(L,0,100,add=TRUE,col="blue")## End(Not run)##################################################### Cox model### piecewise constant hazard, gamma frailty##################################################m <- lvm(y~x+z)rates <- c(0.3,0.5); cuts <- c(0,5)distribution(m,~y+z) <- list(coxExponential.lvm(rate=rates,timecut=cuts),                             loggamma.lvm(rate=1,shape=1))## Not run:     d <- sim(m,2e4,p=c("y~x"=0,"y~z"=0)); d$status <- TRUE    plot(timereg::aalen(survival::Surv(y,status)~x,data=d,                        resample.iid=0,robust=0),spec=1)    L <- approxfun(c(cuts,max(d$y)),f=1,                   cumsum(c(0,rates*diff(c(cuts,max(d$y))))),                   method="linear")    curve(L,0,100,add=TRUE,col="blue")## End(Not run)## Equivalent via transform (here with Aalens additive hazard model)m <- lvm(y~x)distribution(m,~y) <- aalenExponential.lvm(rate=rates,timecut=cuts)distribution(m,~z) <- Gamma.lvm(rate=1,shape=1)transform(m,t~y+z) <- prodsim(m,10)## Shared frailtym <- lvm(c(t1,t2)~x+z)rates <- c(1,0.5); cuts <- c(0,5)distribution(m,~y) <- aalenExponential.lvm(rate=rates,timecut=cuts)distribution(m,~z) <- loggamma.lvm(rate=1,shape=1)## Not run: mets::fast.reshape(sim(m,100),varying="t")## End(Not run)##################################################### General multivariate distributions#################################################### Not run: m <- lvm()distribution(m,~y1+y2,oratio=4) <- VGAM::rbiplackcopksmooth2(sim(m,1e4),rgl=FALSE,theta=-20,phi=25)m <- lvm()distribution(m,~z1+z2,"or1") <- VGAM::rbiplackcopdistribution(m,~y1+y2,"or2") <- VGAM::rbiplackcopsim(m,10,p=c(or1=0.1,or2=4))## End(Not run)m <- lvm()distribution(m,~y1+y2+y3,TRUE) <- function(n,...) rmvn0(n,sigma=diag(3)+1)var(sim(m,100))## Syntax also useful for univariate generators, e.g.m <- lvm(y~x+z)distribution(m,~y,TRUE) <- function(n) rnorm(n,mean=1000)sim(m,5)distribution(m,~y,"m1",0) <- rnormsim(m,5)sim(m,5,p=c(m1=100))##################################################### Regression design in other parameters#################################################### Variance heterogeneitym <- lvm(y~x)distribution(m,~y) <- function(n,mean,x) rnorm(n,mean,exp(x)^.5)if (interactive()) plot(y~x,sim(m,1e3))## Alternaively, calculate the standard error directlyaddvar(m) <- ~sd ## If 'sd' should be part of the resulting data.frameconstrain(m,sd~x) <- function(x) exp(x)^.5distribution(m,~y) <- function(n,mean,sd) rnorm(n,mean,sd)if (interactive()) plot(y~x,sim(m,1e3))## Regression on variance parameterm <- lvm()regression(m) <- y~xregression(m) <- v~x##distribution(m,~v) <- 0 # No stochastic term## Alternative:## regression(m) <- v[NA:0]~xdistribution(m,~y) <- function(n,mean,v) rnorm(n,mean,exp(v)^.5)if (interactive()) plot(y~x,sim(m,1e3))## Regression on shape parameter in Weibull modelm <- lvm()regression(m) <- y ~ z+vregression(m) <- s ~ exp(0.6*x-0.5*z)distribution(m,~x+z) <- binomial.lvm()distribution(m,~cens) <- coxWeibull.lvm(scale=1)distribution(m,~y) <- coxWeibull.lvm(scale=0.1,shape=~s)eventTime(m) <- time ~ min(y=1,cens=0)if (interactive()) {    d <- sim(m,1e3)    require(survival)    (cc <- coxph(Surv(time,status)~v+strata(x,z),data=d))    plot(survfit(cc) ,col=1:4,mark.time=FALSE)}##################################################### Categorical predictor##################################################m <- lvm()## categorical(m,K=3) <- "v"categorical(m,labels=c("A","B","C")) <- "v"regression(m,additive=FALSE) <- y~v## Not run: plot(y~v,sim(m,1000,p=c("y~v:2"=3)))## End(Not run)m <- lvm()categorical(m,labels=c("A","B","C"),p=c(0.5,0.3)) <- "v"regression(m,additive=FALSE,beta=c(0,2,-1)) <- y~v## equivalent to:## regression(m,y~v,additive=FALSE) <- c(0,2,-1)regression(m,additive=FALSE,beta=c(0,4,-1)) <- z~vtable(sim(m,1e4)$v)glm(y~v, data=sim(m,1e4))glm(y~v, data=sim(m,1e4,p=c("y~v:1"=3)))transform(m,v2~v) <- function(x) x=='A'sim(m,10)##################################################### Pre-calculate object##################################################m <- lvm(y~x)m2 <- sim(m,'y~x'=2)sim(m,10,'y~x'=2)sim(m2,10) ## Faster

Spaghetti plot

Description

Spaghetti plot for longitudinal data

Usage

spaghetti(  formula,  data = NULL,  id = "id",  group = NULL,  type = "o",  lty = 1,  pch = NA,  col = 1:10,  alpha = 0.3,  lwd = 1,  level = 0.95,  trend.formula = formula,  tau = NULL,  trend.lty = 1,  trend.join = TRUE,  trend.delta = 0.2,  trend = !is.null(tau),  trend.col = col,  trend.alpha = 0.2,  trend.lwd = 3,  trend.jitter = 0,  legend = NULL,  by = NULL,  xlab = "Time",  ylab = "",  add = FALSE,  ...)

Arguments

Author(s)

Klaus K. Holst

Examples

if (interactive() & requireNamespace("mets")) {K <- 5y <- "y"%++%seq(K)m <- lvm()regression(m,y=y,x=~u) <- 1regression(m,y=y,x=~s) <- seq(K)-1regression(m,y=y,x=~x) <- "b"N <- 50d <- sim(m,N); d$z <- rbinom(N,1,0.5)dd <- mets::fast.reshape(d); dd$num <- dd$num+3spaghetti(y~num,dd,id="id",lty=1,col=Col(1,.4),          trend.formula=~factor(num),trend=TRUE,trend.col="darkblue")dd$num <- dd$num+rnorm(nrow(dd),sd=0.5) ## Unbalancespaghetti(y~num,dd,id="id",lty=1,col=Col(1,.4),          trend=TRUE,trend.col="darkblue")spaghetti(y~num,dd,id="id",lty=1,col=Col(1,.4),           trend.formula=~num+I(num^2),trend=TRUE,trend.col="darkblue")}

Stack estimating equations

Description

Stack estimating equations (two-stage estimator)

Usage

## S3 method for class 'estimate'stack(  x,  model2,  D1u,  inv.D2u,  propensity,  dpropensity,  U,  keep1 = FALSE,  propensity.arg,  estimate.arg,  na.action = na.pass,  ...)

Arguments

Examples

m <- lvm(z0~x)Missing(m, z ~ z0) <- r~xdistribution(m,~x) <- binomial.lvm()p <- c(r=-1,'r~x'=0.5,'z0~x'=2)beta <- p[3]/2d <- sim(m,500,p=p,seed=1)m1 <- estimate(r~x,data=d,family=binomial)d$w <- d$r/predict(m1,type="response")m2 <- estimate(z~1, weights=w, data=d)(e <- stack(m1,m2,propensity=TRUE))

For internal use

Description

For internal use

Author(s)

Klaus K. Holst

Extract subset of latent variable model

Description

Extract measurement models or user-specified subset of model

Usage

## S3 method for class 'lvm'subset(x, vars, ...)

Arguments

Value

Alvm-object.

Author(s)

Klaus K. Holst

Examples

m <- lvm(c(y1,y2)~x1+x2)subset(m,~y1+x1)

Summary method for 'sim' objects

Description

Summary method for 'sim' objects

Usage

## S3 method for class 'sim'summary(  object,  estimate = NULL,  se = NULL,  confint = !is.null(se) && !is.null(true),  true = NULL,  fun,  names = NULL,  unique.names = TRUE,  minimal = FALSE,  level = 0.95,  quantiles = c(0, 0.025, 0.5, 0.975, 1),  ...)

Arguments

Time-dependent parameters

Description

Add time-varying covariate effects to model

Usage

timedep(object, formula, rate, timecut, type = "coxExponential.lvm", ...)

Arguments

Author(s)

Klaus K. Holst

Examples

## Piecewise constant hazardm <- lvm(y~1)m <- timedep(m,y~1,timecut=c(0,5),rate=c(0.5,0.3))## Not run: d <- sim(m,1e4); d$status <- TRUEdd <- mets::lifetable(Surv(y,status)~1,data=d,breaks=c(0,5,10));exp(coef(glm(events ~ offset(log(atrisk)) + -1 + interval, dd, family=poisson)))## End(Not run)## Piecewise constant hazard and time-varying effect of z1m <- lvm(y~1)distribution(m,~z1) <- Binary.lvm(0.5)R <- log(cbind(c(0.2,0.7,0.9),c(0.5,0.3,0.3)))m <- timedep(m,y~z1,timecut=c(0,3,5),rate=R)## Not run: d <- sim(m,1e4); d$status <- TRUEdd <- mets::lifetable(Surv(y,status)~z1,data=d,breaks=c(0,3,5,Inf));exp(coef(glm(events ~ offset(log(atrisk)) + -1 + interval+z1:interval, dd, family=poisson)))## End(Not run)## Explicit simulation of time-varying effectsm <- lvm(y~1)distribution(m,~z1) <- Binary.lvm(0.5)distribution(m,~z2) <- binomial.lvm(p=0.5)#variance(m,~m1+m2) <- 0#regression(m,m1[m1:0] ~ z1) <- log(0.5)#regression(m,m2[m2:0] ~ z1) <- log(0.3)regression(m,m1 ~ z1,variance=0) <- log(0.5)regression(m,m2 ~ z1,variance=0) <- log(0.3)intercept(m,~m1+m2) <- c(-0.5,0)m <- timedep(m,y~m1+m2,timecut=c(0,5))## Not run: d <- sim(m,1e5); d$status <- TRUEdd <- mets::lifetable(Surv(y,status)~z1,data=d,breaks=c(0,5,Inf))exp(coef(glm(events ~ offset(log(atrisk)) + -1 + interval + interval:z1, dd, family=poisson)))## End(Not run)

Converts strings to formula

Description

Converts a vector of predictors and a vector of responses (characters) i#ntoa formula expression.

Usage

toformula(y = ".", x = ".")

Arguments

Value

An object of classformula

Author(s)

Klaus K. Holst

See Also

as.formula,

Examples

toformula(c("age","gender"), "weight")

Trace operator

Description

Calculates the trace of a square matrix.

Usage

tr(x, ...)

Arguments

Value

numeric

Author(s)

Klaus K. Holst

See Also

crossprod,tcrossprod

Examples

tr(diag(1:5))

Trim string of (leading/trailing/all) white spaces

Description

Trim string of (leading/trailing/all) white spaces

Usage

trim(x, all = FALSE, ...)

Arguments

Author(s)

Klaus K. Holst

Twin menarche data

Description

Simulated data

Format

data.frame

Source

Simulated

Two-stage estimator

Description

Generic function.

Usage

twostage(object, ...)

Arguments

See Also

twostage.lvm twostage.lvmfit twostage.lvm.mixture twostage.estimate

Two-stage estimator (non-linear SEM)

Description

Two-stage estimator for non-linear structural equation models

Usage

## S3 method for class 'lvmfit'twostage(  object,  model2,  data = NULL,  predict.fun = NULL,  id1 = NULL,  id2 = NULL,  all = FALSE,  formula = NULL,  std.err = TRUE,  ...)

Arguments

Examples

m <- lvm(c(x1,x2,x3)~f1,f1~z,         c(y1,y2,y3)~f2,f2~f1+z)latent(m) <- ~f1+f2d <- simulate(m,100,p=c("f2,f2"=2,"f1,f1"=0.5),seed=1)## Full MLEee <- estimate(m,d)## Manual two-stage## Not run: m1 <- lvm(c(x1,x2,x3)~f1,f1~z); latent(m1) <- ~f1e1 <- estimate(m1,d)pp1 <- predict(e1,f1~x1+x2+x3)d$u1 <- pp1[,]d$u2 <- pp1[,]^2+attr(pp1,"cond.var")[1]m2 <- lvm(c(y1,y2,y3)~eta,c(y1,eta)~u1+u2+z); latent(m2) <- ~etae2 <- estimate(m2,d)## End(Not run)## Two-stagem1 <- lvm(c(x1,x2,x3)~f1,f1~z); latent(m1) <- ~f1m2 <- lvm(c(y1,y2,y3)~eta,c(y1,eta)~u1+u2+z); latent(m2) <- ~etapred <- function(mu,var,data,...)    cbind("u1"=mu[,1],"u2"=mu[,1]^2+var[1])(mm <- twostage(m1,model2=m2,data=d,predict.fun=pred))if (interactive()) {    pf <- function(p) p["eta"]+p["eta~u1"]*u + p["eta~u2"]*u^2    plot(mm,f=pf,data=data.frame(u=seq(-2,2,length.out=100)),lwd=2)} ## Reduce test timing## Splinesf <- function(x) cos(2*x)+x+-0.25*x^2m <- lvm(x1+x2+x3~eta1, y1+y2+y3~eta2, latent=~eta1+eta2)functional(m, eta2~eta1) <- fd <- sim(m,500,seed=1,latent=TRUE)m1 <- lvm(x1+x2+x3~eta1,latent=~eta1)m2 <- lvm(y1+y2+y3~eta2,latent=~eta2)mm <- twostage(m1,m2,formula=eta2~eta1,type="spline")if (interactive()) plot(mm)nonlinear(m2,type="quadratic") <- eta2~eta1a <- twostage(m1,m2,data=d)if (interactive()) plot(a)kn <- c(-1,0,1)nonlinear(m2,type="spline",knots=kn) <- eta2~eta1a <- twostage(m1,m2,data=d)x <- seq(-3,3,by=0.1)y <- predict(a, newdata=data.frame(eta1=x))if (interactive()) {  plot(eta2~eta1, data=d)  lines(x,y, col="red", lwd=5)  p <- estimate(a,f=function(p) predict(a,p=p,newdata=x))$coefmat  plot(eta2~eta1, data=d)  lines(x,p[,1], col="red", lwd=5)  confband(x,lower=p[,3],upper=p[,4],center=p[,1], polygon=TRUE, col=Col(2,0.2))  l1 <- lm(eta2~splines::ns(eta1,knots=kn),data=d)  p1 <- predict(l1,newdata=data.frame(eta1=x),interval="confidence")  lines(x,p1[,1],col="green",lwd=5)  confband(x,lower=p1[,2],upper=p1[,3],center=p1[,1], polygon=TRUE, col=Col(3,0.2))} ## Reduce test timing## Not run:  ## Reduce timing ## Cross-validation example ma <- lvm(c(x1,x2,x3)~u,latent=~u) ms <- functional(ma, y~u, value=function(x) -.4*x^2) d <- sim(ms,500)#,seed=1) ea <- estimate(ma,d) mb <- lvm() mb1 <- nonlinear(mb,type="linear",y~u) mb2 <- nonlinear(mb,type="quadratic",y~u) mb3 <- nonlinear(mb,type="spline",knots=c(-3,-1,0,1,3),y~u) mb4 <- nonlinear(mb,type="spline",knots=c(-3,-2,-1,0,1,2,3),y~u) ff <- lapply(list(mb1,mb2,mb3,mb4),      function(m) function(data,...) twostage(ma,m,data=data,st.derr=FALSE)) a <- cv(ff,data=d,rep=1) a## End(Not run)

Cross-validated two-stage estimator

Description

Cross-validated two-stage estimator for non-linear SEM

Usage

twostageCV(  model1,  model2,  data,  control1 = list(trace = 0),  control2 = list(trace = 0),  knots.boundary,  nmix = 1:4,  df = 1:9,  fix = TRUE,  std.err = TRUE,  nfolds = 5,  rep = 1,  messages = 0,  ...)

Arguments

Examples

 ## Reduce Ex.Timings##'m1 <- lvm( x1+x2+x3 ~ u, latent= ~u)m2 <- lvm( y ~ 1 )m <- functional(merge(m1,m2), y ~ u, value=function(x) sin(x)+x)distribution(m, ~u1) <- uniform.lvm(-6,6)d <- sim(m,n=500,seed=1)nonlinear(m2) <- y~u1if (requireNamespace('mets', quietly=TRUE)) {  set.seed(1)  val <- twostageCV(m1, m2, data=d, std.err=FALSE, df=2:6, nmix=1:2,                  nfolds=2)  val}

Extract variable names from latent variable model

Description

Extract exogenous variables (predictors), endogenous variables (outcomes),latent variables (random effects), manifest (observed) variables from alvm object.

Usage

vars(x,...)endogenous(x,...)exogenous(x,...)manifest(x,...)latent(x,...)## S3 replacement method for class 'lvm'exogenous(x, xfree = TRUE,...) <- value## S3 method for class 'lvm'exogenous(x,variable,latent=FALSE,index=TRUE,...)## S3 replacement method for class 'lvm'latent(x,clear=FALSE,...) <- value

Arguments

Details

vars returns all variables of thelvm-object includingmanifest and latent variables. Similarilymanifest andlatentreturns the observered resp. latent variables of the model.exogenous returns all manifest variables without parents, e.g.covariates in the model, however the argumentlatent=TRUE can be usedto also include latent variables without parents in the result. Pr. defaultlava will not include the parameters of the exogenous variables inthe optimisation routine during estimation (likelihood of the remainingobservered variables conditional on the covariates), however this behaviourcan be altered via the assignment functionexogenous<- tellinglava which subset of (valid) variables to condition on. Finallylatent returns a vector with the names of the latent variables inx. The assigment functionlatent<- can be used to change thelatent status of variables in the model.

Value

Vector of variable names.

Author(s)

Klaus K. Holst

See Also

endogenous,manifest,latent,exogenous,vars

Examples

g <- lvm(eta1 ~ x1+x2)regression(g) <- c(y1,y2,y3) ~ eta1latent(g) <- ~eta1endogenous(g)exogenous(g)identical(latent(g), setdiff(vars(g),manifest(g)))

vec operator

Description

vec operator

Usage

vec(x, matrix = FALSE, sep = ".", ...)

Arguments

Details

Convert array into vector

Author(s)

Klaus Holst

Wait for user input (keyboard or mouse)

Description

Wait for user input (keyboard or mouse)

Usage

wait()

Author(s)

Klaus K. Holst

Weighted K-means

Description

Weighted K-means via Lloyd's algorithm

Usage

wkm(  x,  mu,  data,  weights = rep(1, NROW(x)),  iter.max = 20,  n.start = 5,  init = "kmpp",  ...)

Arguments

Author(s)

Klaus K. Holst

Wrap vector

Description

Wrap vector

Usage

wrapvec(x, delta = 0L, ...)

Arguments

Examples

wrapvec(5,2)

Regression model for binomial data with unkown group of immortals

Description

Regression model for binomial data with unkown group of immortals (zero-inflated binomial regression)

Usage

zibreg(  formula,  formula.p = ~1,  data,  family = stats::binomial(),  offset = NULL,  start,  var = "hessian",  ...)

Arguments

Author(s)

Klaus K. Holst

Examples

## Simulationn <- 2e3x <- runif(n,0,20)age <- runif(n,10,30)z0 <- rnorm(n,mean=-1+0.05*age)z <- cut(z0,breaks=c(-Inf,-1,0,1,Inf))p0 <- lava:::expit(model.matrix(~z+age) %*% c(-.4, -.4, 0.2, 2, -0.05))y <- (runif(n)<lava:::tigol(-1+0.25*x-0*age))*1u <- runif(n)<p0y[u==0] <- 0d <- data.frame(y=y,x=x,u=u*1,z=z,age=age)head(d)## Estimatione0 <- zibreg(y~x*z,~1+z+age,data=d)e <- zibreg(y~x,~1+z+age,data=d)compare(e,e0)ePD(e0,intercept=c(1,3),slope=c(2,6))B <- rbind(c(1,0,0,0,20),           c(1,1,0,0,20),           c(1,0,1,0,20),           c(1,0,0,1,20))prev <- summary(e,pr.contrast=B)$prevalencex <- seq(0,100,length.out=100)newdata <- expand.grid(x=x,age=20,z=levels(d$z))fit <- predict(e,newdata=newdata)plot(0,0,type="n",xlim=c(0,101),ylim=c(0,1),xlab="x",ylab="Probability(Event)")count <- 0for (i in levels(newdata$z)) {  count <- count+1  lines(x,fit[which(newdata$z==i)],col="darkblue",lty=count)}abline(h=prev[3:4,1],lty=3:4,col="gray")abline(h=prev[3:4,2],lty=3:4,col="lightgray")abline(h=prev[3:4,3],lty=3:4,col="lightgray")legend("topleft",levels(d$z),col="darkblue",lty=seq_len(length(levels(d$z))))