Package 'lava'

Description

For matrices a block-diagonal matrix is created. For all other data types he operator is a wrapper of paste.

Usage

x %++% y

Arguments

Details

Concatenation operator

Author(s)

Klaus K. Holst

See Also

blockdiag, paste, cat,

Examples

## Block diagonal
matrix(rnorm(25),5)%++%matrix(rnorm(25),5)
## String concatenation
"Hello "%++%" World"
## Function composition
f <- log %++% exp
f(2)

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)

Description

Generic method for adding variables to model object

Usage

addvar(x, ...)

Arguments

Author(s)

Klaus K. Holst

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)

Description

Generic method for labeling elements of an object

Usage

baptize(x, ...)

Arguments

Author(s)

Klaus K. Holst

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

Examples

## ---------------------------------------------------------------
## binomial.rd: constant risk-difference model
##   P(Y=1|Z=1) - P(Y=1|Z=0) = tanh(lp)
## ---------------------------------------------------------------
m <- lvm()
regression(m) <- z ~ x
regression(m) <- lp ~ x
regression(m) <- op ~ x
intercept(m, ~lp) <- 0.4   ## constant linear predictor for RD
intercept(m, ~op) <- 0     ## odds product = exp(0) = 1
distribution(m, ~lp) <- normal.lvm(sd = 0)
distribution(m, ~op) <- normal.lvm(sd = 0)
m <- binomial.rd(m, response = "y", exposure = "z",
                 target.model = "lp", nuisance.model = "op")
set.seed(1)
d <- sim(m, n = 2000)
## Empirical risk difference should be close to tanh(0.4)
mean(d$y[d$z == 1]) - mean(d$y[d$z == 0])
tanh(0.4)

## Formula interface: response ~ exposure | target | nuisance
m2 <- lvm()
regression(m2) <- z ~ x
regression(m2) <- lp ~ x
regression(m2) <- op ~ x
m2 <- binomial.rd(m2, y ~ z | lp | op)

## ---------------------------------------------------------------
## binomial.rr: constant relative-risk model
##   log(P(Y=1|Z=1) / P(Y=1|Z=0)) = lp
## ---------------------------------------------------------------
m <- lvm()
regression(m) <- z ~ x
regression(m) <- lp ~ x
regression(m) <- op ~ x
intercept(m, ~lp) <- log(1.5)   ## constant log relative-risk
intercept(m, ~op) <- 0          ## odds product = 1
distribution(m, ~lp) <- normal.lvm(sd = 0)
distribution(m, ~op) <- normal.lvm(sd = 0)
m <- binomial.rr(m, response = "y", exposure = "z",
                 target.model = "lp", nuisance.model = "op")
set.seed(1)
d <- sim(m, n = 2000)
## Empirical log-RR should be close to log(1.5)
log(mean(d$y[d$z == 1]) / mean(d$y[d$z == 0]))
log(1.5)

Description

Combine matrices to block diagonal structure

Usage

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

Arguments

Author(s)

Klaus K. Holst

Examples

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

Description

Bone Mineral Density Data consisting of 112 girls randomized to receive calcium 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

Description

Description

Format

data.frame

Description

Generic method for calculating bootstrap statistics

Usage

bootstrap(x, ...)

Arguments

Author(s)

Klaus K. Holst

See Also

bootstrap.lvm bootstrap.lvmfit

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

A bootstrap.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.Timings
B <- bootstrap(e,R=50,mc.cores=1)
B

Description

Simulated data

Format

data.frame

Source

Simulated

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)

Description

Bone Mineral Density Data consisting of 112 girls randomized to receive calcium 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 at visiting times 'visit', and age in years 'age' at the time of visit. The bone mineral density variable is 'bmd' (g/cm^2).

Source

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

Description

Generic cancel method

Usage

cancel(x, ...)

Arguments

Author(s)

Klaus K. Holst

Description

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

Usage

children(object, ...)

Arguments

Author(s)

Klaus K. Holst

Description

Extension of the identify 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 is terminated by pressing any mouse button other than the first. For the 'quartz' device the process is terminated by pressing either the 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))
}

Description

Given $k$ hypotheses $H_1, \ldots, H_k$ all $2^k-1$ intersection hypotheses are calculated and adjusted p-values are obtained for $H_j$ is calculated as the max $p$-value of all intersection hypotheses containing $H_j$. Example, for $k=3$, the adjusted $p$-value for $H_1$ will be obtained from $\{(H_1, H_2, H_3), (H_1,H_2), (H_1,H_3), (H_1)\}$.

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 with an argument index which is a integer vector specifying which subcomponents of object to test. The ellipsis argument can be any other arguments used in the test function. The function test_wald() is an example of valid test function (which has an additional argument null in reference to the above mentioned ellipsis arguments).

References

Marcus, R; Peritz, E; Gabriel, KR (1976). "On closed testing procedures with special reference to ordered analysis of 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") <- 1
variance(m,endogenous(m)) <- 1
set.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, test = test_wald, null = 0)
adj
adj$p.value
summary(adj)

Description

This function transforms a standard color (e.g. "red") into an transparent 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 in hexadecimal (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)

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)

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)))

Description

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

Usage

commutation(m, n = m)

Arguments

Author(s)

Klaus K. Holst

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) <- ~eta
regression(m) <- eta ~ x
m2 <- 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=0
compare(e,scoretest=y3~x) ## Score-test, H0: y3~x=0
compare(e2,par=c("y3~x")) ## Wald-test, H0: y3~x=0

B <- 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] <- 1
compare(e2,contrast=B)

compare(e,scoretest=list(y3~x,y2~x))

Description

Estimate parameters in a probit latent variable model via a composite likelihood decomposition.

Usage

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

Arguments

Value

An object of class estimate.complik inheriting methods from lvm

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+y3
d <- sim(m,50,seed=1)
if (requireNamespace("mets", quietly=TRUE)) {
   e1 <- complik(m,d,control=list(trace=1),type="all")
}

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,
  alpha = 1,
  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 <- 20
est <- 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),] <- NA
rownames(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="darkred", alpha=0.1)

set.seed(1)
k <- 10
x <- 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 - 1
yu <- y + 1
plot_region(x, y, yl, yu)
plot_region(x, y, yl, yu, type='s', col="darkblue", add=TRUE)

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} \pm z_{\alpha/2}*\hat\sigma_{\hat\theta}$

or profile likelihood confidence limits, defined as the set of value $\tau$:

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

where $q_{\alpha}$ is the $\alpha$ fractile of the $\chi^2_1$ distribution, and $\hat\theta_{\tau}$ are obtained by maximizing the log-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.timings
B <- bootstrap(e,R=50)
B

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 <- 200
x <- seq(0,6,length.out=n)
delta <- 3
ss <- exp(-1+1.5*cos((x-delta)))
ee <- rnorm(n,sd=ss)
y <- (x-delta)+3*cos(x+4.5-delta)+ee
d <- 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))
}

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,...) <- value

constraints(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 associations between covariates and latent or observed variables in the model (non-linear regression).

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 arguments args can be covariates in the model, allowing the specification of non-linear regression models. As an example the non-linear regression model

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

where $\Phi$ denotes the standard normal cumulative distribution function, can be defined as

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

Next we add three new parameters using the parameter assigment function:

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

The intercept of $Y$ is defined as mu

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

And finally the newly added intercept parameter mu is defined as the appropiate non-linear function of $\alpha$, $\nu$ and $\beta$:

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

The constraints function can be used to show the estimated non-linear parameter constraints of an estimated model object (lvmfit or multigroupfit). Calling constrain with no additional arguments beyound x will return a list of the functions and parameter names defining the non-linear restrictions.

The gradient function can optionally be added as an attribute grad to the return value of the function defined by value. In this case the analytical derivatives will be calculated via the chain rule when evaluating the corresponding score function of the log-likelihood. If the gradient attribute is omitted the chain rule will be applied on a numeric approximation of the gradient.

Value

A lvm 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 model
e <- 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 model
logLik(e,p1)
## Likelihood of restricted model (MLE)
logLik(e1)

## Non-linear regression

## Simulate data
m <- lvm(c(y1,y2)~f(x,0)+f(eta,1))
latent(m) <- ~eta
covariance(m,~y1+y2) <- "v"
intercept(m,~y1+y2) <- "mu"
covariance(m,~eta) <- "zeta"
intercept(m,~eta) <- 0
set.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 parameters
constrain(m, mu ~ x + alpha + nu + gamma) <- function(x) x[4]*pnorm(x[3]+x[1]*x[2])
 ## Reduce Ex.Timings
e <- estimate(m,d,control=list(trace=1,constrain=TRUE))
constraints(e,data=d)
## Plot model-fit
plot(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
m1 <- lvm(y ~ f(x,beta)+f(z,beta2))
m2 <- lvm(y ~ f(x,psi) + z)
# simulate data from the two models
d1 <- sim(m1,500)
d2 <- sim(m2,500)
# Add 'non'-linear parameter constraint
constrain(m2,psi ~ beta2) <- function(x) x
## Add parameter beta2 to model 2, now beta2 exists in both models
parameter(m2) <- ~ beta2
ee <- 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

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)

Description

Generic correlation method

Usage

correlation(x, ...)

Arguments

Author(s)

Klaus K. Holst

Description

Define covariances between residual terms in a lvm-object.

Usage

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

Arguments

Details

The covariance function is used to specify correlation structure between 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 a covariance parameter $r = cov(\epsilon_1,\epsilon_2)$, we execute the following code

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

The special function f and its second argument could be omitted thus assigning an unique parameter the covariance between y1 and y2.

Similarily the marginal variance of the two response variables can be fixed to 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 right handside expression of the assigment function covariance<- if the first (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 to NA (see also the regression method).

The function covariance (called without additional arguments) can be used to inspect the covariance constraints of a lvm-object.

Value

A lvm-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 y2
covariance(m) <- y1~y2
covariance(m) <- c(y1,y2)~f(v) ## Same marginal variance
covariance(m) ## Examine covariance structure

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):n
csplit(iris[1:10,],0.5)

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")
}

Description

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

Format

data.frame

Source

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

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

Description

Returns device-coordinates and plot-region

Usage

devcoords()

Value

A list with elements

• dev.x1: device left x-coordinate

• dev.x2: device right x-coordinate

• dev.y1: device bottom y-coordinate

• dev.y2: device top y-coordinate

• fig.x1: plot left x-coordinate

• fig.x2: plot right x-coordinate

• fig.y1: plot bottom y-coordinate

• fig.y2: plot top y-coordinate

Author(s)

Klaus K. Holst

Description

Calculate prevalence, sensitivity, specificity, and positive and negative 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 and test in rows. Data.frame should be with test in the first column 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)

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+v or ~x+y|z+v With everything on the rhs of the bar defining the variables 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)

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

Description

Estimate parameters for the sample mean, variance, and quantiles

Usage

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

Arguments

Description

Primary tool for obtaining parameter estimates with robust (sandwich) standard errors, applying the delta method, and testing linear hypotheses. The function returns an object of class estimate which serves as a general container for parameter estimates and their influence functions (IFs). Three calling conventions are supported:

Usage

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

Arguments

Details

• estimate(x, ...) – extract estimates from a model object

• estimate(coef=, IC=, ...) – construct from coefficients and IF matrix

• estimate(coef=, vcov=, ...) – construct from coefficients and covariance matrix

Influence functions and robust standard errors

An estimator $\widehat{\theta}$ is regular and asymptotically linear (RAL) when it admits the iid decomposition

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

where $\mathrm{IC}$ is the unique influence function satisfying $E\{\mathrm{IC}(Z; P)\} = 0$. By the central limit theorem

$\sqrt{n}(\widehat{\theta}-\theta) \overset{d}{\longrightarrow} N(0,\; \mathrm{Var}\{\mathrm{IC}(Z; P)\})$

and the asymptotic variance is consistently estimated by the empirical variance of the plugin IF estimate, yielding robust (sandwich) standard errors. The estimated IF can be extracted with the IC method.

Parameter transformations (delta method)

When f is a function $\phi: R^p \to R^m$, the delta method is applied:

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

Derivatives are computed numerically via numDeriv::jacobian unless the function returns an attribute "grad" with the analytic Jacobian.

Alternatively, estimate objects support direct arithmetic operations (e.g., a * b, exp(a), a^b) which apply the delta method with exact (analytical) derivatives computed automatically. This influence function calculus allows building complex transformations from simple building blocks without numerical differentiation. See the last example section ("influence function calculus") and vignette("influencefunction", package = "lava") for details.

Averaging and marginalization

When average = TRUE and f(p, data) depends on covariates, the target parameter is the standardized (marginalized) estimate $\Psi = E\{f(X;\theta)\}$. The IF for the averaged estimate accounts for both the empirical averaging and parameter estimation uncertainty:

$\mathrm{IC}_\Psi(Z; P) = f(X;\theta) - \Psi + [E\nabla_\theta f(X;\theta)]\,\phi(Z; P)$

When subset is also specified, the average is conditioned on the subpopulation, yielding a conditional marginalized estimate.

Cluster-robust standard errors

When id is supplied, the per-observation IF contributions are summed within clusters (when stack = TRUE), producing the cluster-level IF $\widetilde{\mathrm{IC}}(Z_i; P) = \sum_{k=1}^{N_i} \frac{n}{N}\mathrm{IC}(Z_{ik}; P)$. The resulting variance estimate is equivalent to the GEE working independence sandwich estimator.

For full theoretical background and worked examples see vignette("influencefunction", package = "lava").

See Also

estimate.array, merge.estimate, contr, parsedesign, pairwise.diff, summary.estimate, coef.estimate, vcov.estimate, transform.estimate, labels.estimate, IC.estimate

Examples

## Simulation from logistic regression model
m <- 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())

## LRT
estimate(g, g0)

## Plain estimates (robust standard errors)
estimate(g)

## Testing contrasts
estimate(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 syntax
estimate(g, "z", "z"-"x", 2*"z"-3*"x")
estimate(g, "?")  ## Wildcards
estimate(g, "*Int*", "z")
estimate(g, "1", "2"-"3", null = c(0,1))
estimate(g, 2, 3)

## Usual (non-robust) confidence intervals
estimate(g, vcov=TRUE)
estimate(g, vcov=vcov(g))

## Transformations
estimate(g, function(p) p[1]+p[2])

## Multiple parameters
e <- estimate(g, function(p) c(p[1]+p[2], p[1]*p[2]))
e
vcov(e)

## Label new parameters
estimate(g, function(p) list("a1"=p[1]+p[2], "b1"=p[1]*p[2]))
#'
## Multiple group
m <- lvm(y~x)
m <- baptize(m)
d2 <- d1 <- sim(m,50,seed=1)
e <- estimate(list(m,m),list(d1,d2))
estimate(e) ## Wrong
ee <- estimate(e, id=rep(seq(nrow(d1)), 2)) ## Clustered
ee
estimate(lm(y~x,d1))

## Marginalize
f <- function(p,data)
  list(p0=expit(p["(Intercept)"] + p["z"]*data[,"z"]),
       p1=expit(p["(Intercept)"] + p["x"] + p["z"]*data[,"z"]))
e <- estimate(g, f, average=TRUE)
e
estimate(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=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-numbers
id1 <- 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. decomposition
IC(estimate(l1,id=id1,stack=FALSE))
IC(estimate(l1,id=id1))

## Combined i.i.d. decomposition
e1 <- 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


# ------ influence function calculus -------
a <- estimate(coef = c("a" = 0.5), IC = scale(rnorm(10), scale=FALSE), id = 1:10)
b <- estimate(coef = c("b" = 0.8), IC = scale(rnorm(10), scale=FALSE), id = 1:10)

e <- c(a, b) # merge
merge(a, b)
c(e1=a, b) # naming of par
labels(e, c("p1", "p2")) # renaming parameters
e["a"] # subset
subset(e, "a")

# pipes
# c(a, b) |>
#  transform(function(x) x^2) |>
#  subset("a") |>
#  labels("sq")

# Parameter transformation with automatic calculation of derivatives
a * b
(3 * cos(a) / sqrt(b) + 1) / a
expit(c(a,b))
c(sum=sum(e), sum2=a+b,
  prod=prod(e), prod2=a*b)
e %*% e # inner prod.
c(1, 2) %*% e
c(pow = a^b)
a^c(0.5, 2)
c(b=e["a"] * e["b"] / a, also.b=e["b"])

B <- rbind(c(1,-1), c(1,0), c(0,1))
B %*% e
e == 1 # wald-test, null-hypothesis H0: b=1
e == c(1,2)
B %*% e == 1

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 procedures is parsed via the control argument. By default Maximum Likelihood is used assuming multivariate normal distributed measurement errors. A list with one or more of the following elements is expected:

start:

Starting value. The order of the parameters can be shown by calling coef (with mean=TRUE) on the lvm-object or with plot(..., labels=TRUE). Note that this requires a check that it is actual the model being estimated, as estimate might add additional restriction to the model, e.g. through the fix and exo.fix arguments. The lvm-object of a fitted model can be extracted with the Model-function.

starterfun:

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

estimator:

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

meanstructure

Logical variable indicating whether to fit model with meanstructure.

method:

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

control:

Parameters passed to the optimizer (default stats::nlminb).

tol:

Tolerance of optimization constraints on lower limit of variance parameters.

Value

A lvmfit-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 example
m <- lvm(list(y~v1+v2+v3+v4,c(v1,v2,v3,v4)~x))
covariance(m) <- v1~v2+v3+v4
dd <- sim(m,10000) ## Simulate 10000 observations from model
e <- estimate(m, dd) ## Estimate parameters
e

## Using just sufficient statistics
n <- nrow(dd)
e0 <- estimate(m,data=list(S=cov(dd)*(n-1)/n,mu=colMeans(dd),n=n))
rm(dd)

## Multiple group analysis
m <- lvm()
regression(m) <- c(y1,y2,y3)~u
regression(m) <- u~x
d1 <- 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) <- NA
covariance(mm,~u) <- NA
intercept(mm,~u) <- NA
ee <- estimate(list(mm,mm),list(d1,d2))

## Missing data
d0 <- makemissing(d1,cols=1:2)
e0 <- estimate(m,d0,missing=TRUE)
e0

Description

For example, if the model 'm' includes latent event time variables are 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 model one 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 covariates
m0 <- 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.4
regression(m,from="sbp",to="eventtime") <- -0.01
sim(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.7
regression(ma,from="sbp",to="eventtime") <- -0.008
set.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 data

mc <- 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)