Package 'lavaSearch2'

Description

Generic interface to add links to lvm objects.

Usage

addLink(object, ...)

## S3 method for class 'lvm'
addLink(
  object,
  var1,
  var2,
  covariance,
  all.vars = lava::vars(object),
  warnings = FALSE,
  ...
)

## S3 method for class 'lvm.reduced'
addLink(object, ...)

Arguments

Details

The argument all.vars is useful for lvm.reduce object where the command vars(object) does not return all variables. The command vars(object, xlp = TRUE) must be used instead.

Arguments var1 and var2 are passed to initVarlink.

Examples

library(lava)
set.seed(10)

m <- lvm()
regression(m) <- c(y1,y2,y3)~u
regression(m) <- u~x1+x2
latent(m) <- ~u
m2 <- m

addLink(m, x1 ~ y1, covariance = FALSE)
addLink(m, y1 ~ x1, covariance = FALSE)
coef(addLink(m, y1 ~ y2, covariance = TRUE))

addLink(m2, "x1", "y1", covariance = FALSE)
addLink(m2, "y1", "x1", covariance = FALSE)
newM <- addLink(m, "y1", "y2", covariance = TRUE)
coef(newM)

Description

Graphical display of the bias or type 1 error for the output of calibrateType1.

Usage

## S3 method for class 'calibrateType1'
autoplot(
  object,
  type = "bias",
  plot = TRUE,
  color.threshold = "red",
  type.bias = "absolute",
  alpha = 0.05,
  nrow.legend = NULL,
  name2label = NULL,
  color = NULL,
  keep.method = NULL,
  ...
)

Arguments

Details

Method names:

• p.Ztest

• p.Satt

• p.KR

• p.robustZtest

• p.robustSatt

• p.robustKR

Value

An list containing:

• plot: a ggplot object.

• data: the dataset used to generate the ggplot object.

Description

2D-display of the domain used to compute the integral.

Usage

## S3 method for class 'intDensTri'
autoplot(object, coord.plot = c("x", "y1"), plot = TRUE, ...)

Arguments

Value

A ggplot object.

See Also

intDensTri

Description

Display the value of a coefficient across the steps.

Usage

## S3 method for class 'modelsearch2'
autoplot(
  object,
  param,
  ci = TRUE,
  step = 0:nStep(object),
  conf.level = 0.95,
  plot = TRUE,
  add.0 = TRUE,
  ...
)

Arguments

Value

A list containing

• plot: a ggplot object.

• data: the data used to generate the ggplot object.

Examples

## Not run: 
mSim <- lvm(Y~G+X1+X2+X3+X4+X5)
addvar(mSim) <- ~Z1+Z2

set.seed(10)
df.data <- lava::sim(mSim, 1e2)

mBase <- lvm(Y~G)
addvar(mBase) <- ~X1+X2+X3+X4+X5+Z1+Z2
e.lvm <- estimate(mBase, data = df.data)
res <- modelsearch2(e.lvm, method.p.adjust = "holm", alpha = 0.05)
autoplot(res, param = "Y~G")
autoplot(res, param = c("Y","Y~G"))

## End(Not run)

Description

Adjust the p.values using the quantiles of the max statistic.

Usage

calcDistMaxIntegral(
  statistic,
  iid,
  df,
  iid.previous = NULL,
  quantile.previous = NULL,
  quantile.compute = lava.options()$search.calc.quantile.int,
  alpha,
  cpus = 1,
  cl = NULL,
  trace
)

calcDistMaxBootstrap(
  statistic,
  iid,
  iid.previous = NULL,
  quantile.previous = NULL,
  method,
  alpha,
  cpus = 1,
  cl = NULL,
  n.sim,
  trace,
  n.repmax = 100
)

Arguments

Value

A list containing

• p.adjust: the adjusted p-values.

• z: the rejection threshold.

• Sigma: the correlation matrix between the test statistic.

• correctedLevel: the alpha level corrected for conditioning on previous tests.

Examples

library(mvtnorm)

set.seed(10)
n <- 100
p <- 4
link <- letters[1:p]
n.sim <- 1e3 # number of bootstrap simulations 

#### test - not conditional ####
X.iid <- rmvnorm(n, mean = rep(0,p), sigma = diag(1,p))
colnames(X.iid) <- link
statistic <- setNames(1:p,link)


r1 <- calcDistMaxIntegral(statistic = statistic, iid = X.iid, 
            trace = FALSE, alpha = 0.05, df = 1e6) 

r3 <- calcDistMaxBootstrap(statistic = statistic, iid = X.iid,
            method = "residual",
            trace = FALSE, alpha = 0.05, n.sim = n.sim)

r4 <- calcDistMaxBootstrap(statistic = statistic, iid = X.iid,
            method = "wild",
            trace = FALSE, alpha = 0.05, n.sim = n.sim)

rbind(integration = c(r1$p.adjust, quantile = r1$z),
      bootResidual = c(r3$p.adjust, quantile = r3$z),
      bootWild    = c(r4$p.adjust, quantile = r4$z))

#### test - conditional ####
## Not run: 
Z.iid <- rmvnorm(n, mean = rep(0,p+1), sigma = diag(1,p+1))
seqQuantile <- qmvnorm(p = 0.95, delta = rep(0,p+1), sigma = diag(1,p+1), 
                    tail = "both.tails")$quantile

r1c <- calcDistMaxIntegral(statistic = statistic, iid = X.iid,
            iid.previous = Z.iid, quantile.previous =  seqQuantile, 
            trace = FALSE, alpha = 0.05, df = NULL)

r3c <- calcDistMaxBootstrap(statistic = statistic, iid = X.iid,
            iid.previous = Z.iid, quantile.previous =  seqQuantile, method = "residual",
            trace = FALSE, alpha = 0.05, n.sim = n.sim)

r4c <- calcDistMaxBootstrap(statistic = statistic, iid = X.iid,
            iid.previous = Z.iid, quantile.previous =  seqQuantile, method = "wild",
            trace = FALSE, alpha = 0.05, n.sim = n.sim)

rbind(integration = c(r1c$p.adjust, quantile = r1c$z),
      bootResidual = c(r3c$p.adjust, quantile = r3c$z),
      bootWild    = c(r4c$p.adjust, quantile = r4c$z))

## End(Not run)

Description

Compute the type 1 error after selection [EXPERIMENTAL].

Usage

calcType1postSelection(
  level,
  mu,
  Sigma,
  quantile.previous,
  distribution,
  df,
  n = 10,
  correct = TRUE,
  ...
)

Arguments

Details

The number of tests at the current step (i.e. after selection) is assumed to be one less than the number of tests at the previous step (i.e. before selection).

Arguments mu and Sigma must contain the moments for the vector of test statistics before and after selection (in that order).

Value

[numeric] the type 1 error.

Author(s)

Brice Ozenne

Examples

library(mvtnorm)
n <- 350

#### only 2 tests
Sigma <- rbind(c(1,0,0),c(0,1,1),c(0,1,1))
z2 <- qmvnorm(0.95, mean = rep(0,2), sigma = Sigma[1:2,1:2], tail = "both.tails")$quantile

## no selection since strong effect
mu <- c(10,0,0)
calcType1postSelection(0.95, quantile.previous = z2, distribution = "gaussian",
                        mu = mu, Sigma = Sigma, correct = TRUE)

## strong selection
## Not run: 
mu <- c(0,0,0)
levelC <- calcType1postSelection(0.95, quantile.previous = z2, distribution = "gaussian",
                        mu = mu, Sigma = Sigma)
print(levelC) # more liberal than without selection
calcType1postSelection(levelC, quantile.previous = z2, distribution = "gaussian",
                        mu = mu, Sigma = Sigma, correct = FALSE)

## End(Not run)

#### 3 tests
Sigma <- diag(1,5,5)
Sigma[4,2] <- 1
Sigma[2,4] <- 1
Sigma[5,3] <- 1
Sigma[3,5] <- 1

z2 <- qmvnorm(0.95, mean = mu[1:3], sigma = Sigma[1:3,1:3], tails = "both.tails")$quantile

## no selection since strong effect
## Not run: 
mu <- c(10,0,0,0,0)
calcType1postSelection(0.95, quantile.previous = z2, distribution = "gaussian",
                        mu = mu, Sigma = Sigma, correct = TRUE)

## strong selection
mu <- c(0,0,0,0,0)
levelC <- calcType1postSelection(0.95, quantile.previous = z2,
                        mu = mu, Sigma = Sigma, distribution = "gaussian")
calcType1postSelection(levelC, quantile.previous = z2, distribution = "gaussian",
                        mu = mu, Sigma = Sigma, correct = FALSE)

## End(Not run)

Description

Perform a simulation study over one or several sample size to assess the bias of the estimate and the type 1 error of the Wald test and robust Wald test

Usage

calibrateType1(object, param, n.rep, ...)

## S3 method for class 'lvm'
calibrateType1(
  object,
  param,
  n.rep,
  n,
  correction = TRUE,
  warmup = NULL,
  null = NULL,
  F.test = FALSE,
  cluster = NULL,
  generative.object = NULL,
  generative.coef = NULL,
  true.coef = NULL,
  n.true = 1e+06,
  round.true = 2,
  bootstrap = FALSE,
  n.bootstrap = 1000,
  checkType1 = FALSE,
  checkType2 = FALSE,
  dir.save = NULL,
  label.file = NULL,
  seed = NULL,
  cpus = 1,
  trace = 2,
  ...
)

## S3 method for class 'lvmfit'
calibrateType1(
  object,
  param,
  n.rep,
  correction = TRUE,
  F.test = FALSE,
  bootstrap = FALSE,
  n.bootstrap = 1000,
  seed = NULL,
  trace = 2,
  cpus = 1,
  ...
)

Arguments

Value

An object of class calibrateType1.

Author(s)

Brice Ozenne

See Also

link{autoplot.calibrateType1} for a graphical display of the bias or of the type 1 error.

Examples

## Not run: 
#### simulate data ####
m.Sim <- lvm(c(Y1[mu1:sigma]~1*eta,
               Y2[mu2:sigma]~1*eta,
               Y3[mu3:sigma]~1*eta,
               eta~beta1*Group+beta2*Gender))
latent(m.Sim) <- ~eta
categorical(m.Sim, labels = c("M","F")) <- ~Gender

d <- lava::sim(m.Sim, 1e2)

#### calibrate type 1 error on the estimated model ####
m <- lvm(Y1~eta,
         Y2~eta,
         Y3~eta,
         eta~Group+Gender)
e <- lava::estimate(m, data = d)
res <- calibrateType1(e, param = "eta~Group", n.rep = 100)
res <- calibrateType1(e, param = c("eta~Group","Y1~eta"), F.test = TRUE, n.rep = 100)
res <- calibrateType1(e, param = "eta~Group", n.rep = 100, cpus = 4)
summary(res)

## End(Not run)

Description

Check whether the dataset can be used to fit the lvm object.

Usage

checkData(object, data, trace)

## S3 method for class 'lvm'
checkData(object, data, trace = TRUE)

Arguments

Value

Invisible TRUE or FALSE.

Examples

m <- lvm()
regression(m) <- c(y1,y2,y3)~u
regression(m) <- u~x
latent(m) <- ~u

d <- lava::sim(m,1e2)

try(checkData(m, data = d)) # return an error

checkData(m, data = d[,-4])

try(checkData(m, data = d[,-(3:4)])) # return an error

Description

Remove variables with no link.

Usage

clean(x, ...)

## S3 method for class 'lvm'
clean(x, rm.exo = TRUE, rm.endo = TRUE, rm.latent = TRUE, ...)

Arguments

Examples

m <- lvm()
m <- regression(m, x=paste0("x",1:5),y="y1")
m <- regression(m, x=paste0("x",1:5),y="y2")
covariance(m) <- y1~y2

cancel(m) <- y1 ~ x1
cancel(m) <- y2 ~ x1
clean(m)

m <- lvm(y1 ~ eta + x1, y2 ~ eta, y3 ~ eta + x2)
latent(m) <- ~eta
clean(m)
m
cancel(m) <- y1 ~ eta
cancel(m) <- y2 ~ eta
cancel(m) <- y3 ~ eta
clean(m)

Description

Extract the coefficients from a latent variable model. Similar to lava::compare but with small sample correction.

Usage

coef2(object, as.lava, ...)

## S3 method for class 'lvmfit'
coef2(object, as.lava = TRUE, ssc = lava.options()$ssc, ...)

Arguments

Details

When argument object is a lvmfit object, the method first calls estimate2 and then extract the model coefficients.

Value

A numeric vector named with the names of the coefficients.

See Also

estimate2 to obtain lvmfit2 objects.

Examples

#### simulate data ####
set.seed(10)
dW <- sampleRepeated(10, format = "wide")
set.seed(10)
dL <- sampleRepeated(10, format = "long")
dL$time2 <- paste0("visit",dL$time)

#### latent variable models ####
e.lvm <- estimate(lvm(c(Y1,Y2,Y3) ~ 1*eta + X1, eta ~ Z1), data = dW)
coef(e.lvm)
coef2(e.lvm)
coef2(e.lvm, as.lava = FALSE)

Description

Extract specific types of coefficient from a lvm object: covariance coefficient(s) (coefCov), extra parameter(s) (coefExtra), position in the list of models for each coefficient (coefIndexModel), intercept coefficient(s) (coefIntercept), coefficient(s) that are used as reference (coefRef), regression coefficient(s) (coefReg), variance coefficient(s) (coefVar).

Usage

coefCov(object, value, keep.var, ...)

## S3 method for class 'lvm'
coefCov(object, value = FALSE, keep.var = FALSE, ...)

## S3 method for class 'lvmfit'
coefCov(object, value = FALSE, keep.var = FALSE, ...)

## S3 method for class 'multigroup'
coefCov(object, value = FALSE, keep.var = FALSE, ...)

coefExtra(object, value, ...)

## S3 method for class 'lvm'
coefExtra(object, value = FALSE, ...)

## S3 method for class 'lvmfit'
coefExtra(object, value = FALSE, ...)

## S3 method for class 'multigroup'
coefExtra(object, value = FALSE, ...)

coefIndexModel(object, ...)

## S3 method for class 'lvm'
coefIndexModel(object, ...)

## S3 method for class 'lvmfit'
coefIndexModel(object, ...)

## S3 method for class 'multigroup'
coefIndexModel(object, ...)

## S3 method for class 'multigroupfit'
coefIndexModel(object, ...)

coefIntercept(object, value, ...)

## S3 method for class 'lvm'
coefIntercept(object, value = FALSE, ...)

## S3 method for class 'lvmfit'
coefIntercept(object, value = FALSE, ...)

## S3 method for class 'multigroup'
coefIntercept(object, value = FALSE, ...)

coefRef(object, value, ...)

## S3 method for class 'lvmfit'
coefRef(object, value = FALSE, ...)

coefReg(object, value, ...)

## S3 method for class 'lvm'
coefReg(object, value = FALSE, ...)

## S3 method for class 'lvmfit'
coefReg(object, value = FALSE, ...)

## S3 method for class 'multigroup'
coefReg(object, value = FALSE, ...)

coefVar(object, value, ...)

## S3 method for class 'lvm'
coefVar(object, value = FALSE, ...)

## S3 method for class 'lvmfit'
coefVar(object, value = FALSE, ...)

## S3 method for class 'multigroup'
coefVar(object, value = FALSE, ...)

Arguments

Value

A vector containing the names of the positions of the coefficients.

Examples

#### regression ####
m <- lvm(Y~X1+X2)
e <- estimate(m, lava::sim(m, 1e2))

coefCov(m)
coefCov(m, value = TRUE)

coefCov(m, keep.var = TRUE)
coefCov(m, value = TRUE, keep.var = TRUE)

coefIndexModel(m)
coefIndexModel(e)

coefIntercept(m)
coefIntercept(m, value = TRUE)

coefReg(m)
coefReg(m, value = TRUE)

#### LVM ####
m <- lvm()
regression(m) <- c(y1,y2,y3)~u
regression(m) <- u~x1+x2
latent(m) <- ~u
covariance(m) <- y1~y2

m.Sim <- m
categorical(m.Sim, labels = c("a","b","c")) <- ~x2
e <- estimate(m, lava::sim(m.Sim, 1e2))

coefCov(m)
coefCov(m, value = TRUE) 

coefCov(m, keep.var = TRUE)
coefCov(m, value = TRUE, keep.var = TRUE)

coefExtra(m)

coefIndexModel(m)
coefIndexModel(e)

## additional categorical variable 
categorical(m, labels = as.character(1:3)) <- "X1"

coefExtra(m)
coefExtra(m, value = TRUE)

## additional categorical variable
categorical(m, labels = as.character(1:3)) <- "x1"

coefIntercept(m)
coefIntercept(m, value = TRUE)
coefIntercept(e)

coefReg(e, value = TRUE)

#### multigroup ####
m <- lvm(Y~X1+X2)
eG <- estimate(list(m,m), list(lava::sim(m, 1e2), lava::sim(m, 1e2)))

coefIndexModel(eG)

Description

Extract the type of each coefficient of a lvm object.

Usage

coefType(object, as.lava, ...)

## S3 method for class 'lvm'
coefType(object, as.lava = TRUE, data = NULL, ...)

## S3 method for class 'lvmfit'
coefType(object, as.lava = TRUE, ...)

## S3 method for class 'multigroup'
coefType(object, as.lava = TRUE, ...)

Arguments

Details

A lvm can be written as a measurement model:

$Y_i = \nu + \Lambda \eta_i + K X_i + \epsilon_i$

and a structural model:

$\eta_i = \alpha + B \eta_i + \Gamma X_i + \zeta_i$

where $\Psi$ is the variance covariance matrix of the residuals $\zeta$
and $\Sigma$ is the variance covariance matrix of the residuals $\epsilon$.

coefType either returns the Latin/Greek letter corresponding to the coefficients or it groups them:

• intercept: $\nu$ and $\alpha$.

• regression: $\Lambda$, $K$, $B$, and $\Gamma$.

• covariance: extra-diagonal terms of $\Sigma$ and $\Psi$.

• variance: diagonal of $\Sigma$ and $\Psi$.

A link denotes a relationship between two variables. The coefficient are used to represent the strength of the association between two variable, i.e. the strength of a link. A coefficient may corresponds to the strength of one or several link.

Value

coefType returns a data.frame when as.lava=FALSE:

• name: name of the link

• Y: outcome variable

• X: regression variable in the design matrix (could be a transformation of the original variables, e.g. dichotomization).

• data: original variable

• type: type of link

• value: if TRUE, the value of the link is set and not estimated.

• marginal: if TRUE, the value of the link does not impact the estimation.

• detail: a more detailed description of the type of link (see the details section)

• lava: name of the coefficient in lava

When as.lava=TRUE, coefType returns a named vector containing the type of each coefficient.

Examples

#### regression ####
m <- lvm(Y~X1+X2)
e <- estimate(m, lava::sim(m, 1e2))

coefType(m)
coefType(e)

#### LVM ####
m <- lvm()
regression(m) <- c(y1,y2,y3)~u
regression(m) <- u~x1+x2
latent(m) <- ~u
covariance(m) <- y1~y2

m.Sim <- m
categorical(m.Sim, labels = c("a","b","c")) <- ~x2
e <- estimate(m, lava::sim(m.Sim, 1e2))

coefType(m)
coefType(e)

## additional categorical variables 
categorical(m, labels = as.character(1:3)) <- "X1"

coefType(m, as.lava = FALSE)

#### LVM with constrains ####
m <- lvm(c(Y1~0+1*eta1,Y2~0+1*eta1,Y3~0+1*eta1,
          Z1~0+1*eta2,Z2~0+1*eta2,Z3~0+1*eta2))
latent(m) <- ~eta1 + eta2
e <- estimate(m, lava::sim(m,1e2))

coefType(m)
coefType(e)

#### multigroup ####
m <- lvm(Y~X1+X2)
eG <- estimate(list(m,m), list(lava::sim(m, 1e2), lava::sim(m, 1e2)))
coefType(eG)

Description

Combine formula by outcome

Usage

combineFormula(ls.formula, as.formula = TRUE, as.unique = FALSE)

Arguments

Examples

combineFormula(list(Y~X1,Y~X3+X5,Y1~X2))
lava.options(symbols = c("~",","))
combineFormula(list("Y~X1","Y~X3+X5","Y1~X2"))
lava.options(symbols = c("<-","<->"))
combineFormula(list("Y<-X1","Y<-X3+X5","Y1<-X2"))

combineFormula(list(Y~X1,Y~X3+X1,Y1~X2))
combineFormula(list(Y~X1,Y~X3+X1,Y1~X2), as.formula = FALSE)
combineFormula(list(Y~X1,Y~X3+X1,Y1~X2), as.unique = TRUE)

lava.options(symbols = c("~","~~"))
combineFormula(list("Y~X1","Y~X3","Y1~X2"))

Description

Test Linear Hypotheses using Wald statistics in a latent variable model. Similar to lava::compare but with small sample correction.

Usage

compare2(
  object,
  linfct,
  rhs,
  robust,
  cluster,
  as.lava,
  F.test,
  conf.level,
  ...
)

## S3 method for class 'lvmfit'
compare2(
  object,
  linfct = NULL,
  rhs = NULL,
  robust = FALSE,
  cluster = NULL,
  as.lava = TRUE,
  F.test = TRUE,
  conf.level = 0.95,
  ssc = lava.options()$ssc,
  df = lava.options()$df,
  ...
)

## S3 method for class 'lvmfit2'
compare2(
  object,
  linfct = NULL,
  rhs = NULL,
  robust = FALSE,
  cluster = NULL,
  as.lava = TRUE,
  F.test = TRUE,
  conf.level = 0.95,
  ...
)

## S3 method for class 'lvmfit2'
compare(
  object,
  linfct = NULL,
  rhs = NULL,
  robust = FALSE,
  cluster = NULL,
  as.lava = TRUE,
  F.test = TRUE,
  conf.level = 0.95,
  ...
)

Arguments

Details

The linfct argument and rhs specify the set of linear hypotheses to be tested. They can be written:

$linfct * \theta = rhs$

where $\theta$ is the vector of the model coefficients.
The par argument must contain expression(s) involving the model coefficients. For example "beta = 0" or c("-5*beta + alpha = 3","-alpha") are valid expressions if alpha and beta belong to the set of model coefficients. A contrast matrix and the right hand side will be generated inside the function.

When directly specified, the contrast matrix must contain as many columns as there are coefficients in the model (mean and variance coefficients). Each hypothesis correspond to a row in the contrast matrix.

The rhs vector should contain as many elements as there are row in the contrast matrix.

Value

If as.lava=TRUE an object of class htest. Otherwise a data.frame object.

See Also

createContrast to create contrast matrices.
estimate2 to obtain lvmfit2 objects.

Examples

#### simulate data ####
set.seed(10)
mSim <- lvm(Y~0.1*X1+0.2*X2)
categorical(mSim, labels = c("a","b","c")) <- ~X1
transform(mSim, Id~Y) <- function(x){1:NROW(x)}
df.data <- lava::sim(mSim, 1e2)

#### with lvm ####
m <- lvm(Y~X1+X2)
e.lvm <- estimate(m, df.data)

compare2(e.lvm, linfct = c("Y~X1b","Y~X1c","Y~X2"))
compare2(e.lvm, linfct = c("Y~X1b","Y~X1c","Y~X2"), robust = TRUE)

Description

Extract confidence intervals of the coefficients from a latent variable model. Similar to lava::confint but with small sample correction.

Extract estimate, standard error, confidence intervals and p-values associated to each coefficient of a latent variable model. Similar to lava::confint but with small sample correction.

Usage

confint2(object, robust, cluster, transform, as.lava, conf.level, ...)

## S3 method for class 'lvmfit'
confint2(
  object,
  robust = FALSE,
  cluster = NULL,
  transform = NULL,
  as.lava = TRUE,
  conf.level = 0.95,
  ssc = lava.options()$ssc,
  df = lava.options()$df,
  ...
)

model.tables2(object, robust, cluster, transform, as.lava, conf.level, ...)

Arguments

Details

When argument object is a lvmfit object, the method first calls estimate2 and then extract the confidence intervals.

When argument object is a lvmfit object, the method first calls estimate2 and then extract the confidence intervals.

Value

A data.frame with a row per coefficient.

A data.frame with a row per coefficient.

Examples

#### simulate data ####
set.seed(10)
dW <- sampleRepeated(10, format = "wide")
set.seed(10)
dL <- sampleRepeated(10, format = "long")
dL$time2 <- paste0("visit",dL$time)

#### latent variable models ####
e.lvm <- estimate(lvm(c(Y1,Y2,Y3) ~ 1*eta + X1, eta ~ Z1), data = dW)
confint(e.lvm)
confint2(e.lvm)
confint2(e.lvm, as.lava = FALSE)

Description

Conversion of formula into character string or vice versa

Usage

formula2character(f, type = "formula")

Arguments

Examples

formula2character(Y1~X1+X2)
formula2character(Y1~X1+X2, type = "symbols")

Description

Returns a contrast matrix corresponding an object. The contrast matrix will contains the hypotheses in rows and the model coefficients in columns.

Usage

createContrast(object, ...)

## S3 method for class 'character'
createContrast(object, ...)

## S3 method for class 'lvmfit'
createContrast(object, linfct, ...)

## S3 method for class 'lvmfit2'
createContrast(object, linfct, ...)

## S3 method for class 'list'
createContrast(object, linfct = NULL, ...)

## S3 method for class 'mmm'
createContrast(object, linfct = NULL, ...)

Arguments

Details

One can initialize an empty contrast matrix setting the argumentlinfct to character(0).

Value

A list containing

• contrast [matrix] a contrast matrix corresponding to the left hand side of the linear hypotheses.

• null [vector] the right hand side of the linear hypotheses.

• Q [integer] the rank of the contrast matrix.

• ls.contrast [list, optional] the contrast matrix corresponding to each submodel. Only present when the argument object is a list of models.

Examples

## Simulate data
mSim <- lvm(X ~ Age + Treatment,
            Y ~ Gender + Treatment,
            c(Z1,Z2,Z3) ~ eta, eta ~ treatment,
            Age[40:5]~1)
latent(mSim) <- ~eta
categorical(mSim, labels = c("placebo","SSRI")) <- ~Treatment
categorical(mSim, labels = c("male","female")) <- ~Gender
n <- 1e2
set.seed(10)
df.data <- lava::sim(mSim,n)

## Estimate separate models
lmX <- lava::estimate(lvm(X ~ -1 + Age + Treatment), data = df.data)
lmY <- lava::estimate(lvm(Y ~ -1 + Gender + Treatment), data = df.data)
lvmZ <- lava::estimate(lvm(c(Z1,Z2,Z3) ~ -1 + 1*eta, eta ~ -1 + Treatment), 
                 data = df.data)

## Contrast matrix for a given model
createContrast(lmX, linfct = "X~Age")
createContrast(lmX, linfct = c("X~Age=0","X~Age+5*X~TreatmentSSRI=0"))
createContrast(lmX, linfct = c("X~Age=0","X~Age+5*X~TreatmentSSRI=0"), sep = NULL)
createContrast(lmX, linfct = character(0))

## Contrast matrix for the join model
ls.lvm <- list(X = lmX, Y = lmY, Z = lvmZ)
createContrast(ls.lvm, linfct = "TreatmentSSRI=0")
createContrast(ls.lvm, linfct = "TreatmentSSRI=0", rowname.rhs = FALSE)
createContrast(ls.lvm, linfct = character(0))

## Contrast for multigroup models
m <- lava::lvm(Y~Age+Treatment)
e <- lava::estimate(list(m,m), data = split(df.data, df.data$Gender))
print(coef(e))
createContrast(e, linfct = "Y~TreatmentSSRI@1 - Y~TreatmentSSRI@2 = 0")
createContrast(e, linfct = "Y~TreatmentSSRI@2 - Y~TreatmentSSRI@1 = 0")

Description

Computation of the degrees of freedom of the chi-squared distribution relative to the model-based variance

Usage

dfSigma(contrast, score, vcov, rvcov, dVcov, dRvcov, keep.param, type)

Arguments

References

Wei Pan and Melanie M. Wall, Small-sample adjustments in using the sandwich variance estiamtor in generalized estimating equations. Statistics in medicine (2002) 21:1429-1441.

Description

Test whether a path in the latent variable model correspond to a null effect. Similar to lava::effects but with small sample correction (if any). So far it only work for a single path related two variable composed of one or two edges.

Usage

effects2(object, linfct, robust, cluster, conf.level, ...)

## S3 method for class 'lvmfit'
effects2(
  object,
  linfct,
  robust = FALSE,
  cluster = NULL,
  conf.level = 0.95,
  to = NULL,
  from = NULL,
  df = lava.options()$df,
  ssc = lava.options()$ssc,
  ...
)

## S3 method for class 'lvmfit2'
effects2(
  object,
  linfct,
  robust = FALSE,
  cluster = NULL,
  conf.level = 0.95,
  to = NULL,
  from = NULL,
  ...
)

## S3 method for class 'lvmfit2'
effects(
  object,
  linfct,
  robust = FALSE,
  cluster = NULL,
  conf.level = 0.95,
  to = NULL,
  from = NULL,
  ...
)

Arguments

Details

When argument object is a lvmfit object, the method first calls estimate2 and then extract the confidence intervals.

Value

A data.frame with a row per path.

Description

Extract data from a latent variable model.

Usage

extractData(object, design.matrix, as.data.frame, envir, rm.na)

## S3 method for class 'lvmfit'
extractData(
  object,
  design.matrix = FALSE,
  as.data.frame = TRUE,
  envir = environment(),
  rm.na = TRUE
)

Arguments

Value

a dataset.

Examples

#### simulate data ####
set.seed(10)
n <- 101

Y1 <- rnorm(n, mean = 0)
Y2 <- rnorm(n, mean = 0.3)
Id <- findInterval(runif(n), seq(0.1,1,0.1))
data.df <- rbind(data.frame(Y=Y1,G="1",Id = Id),
           data.frame(Y=Y2,G="2",Id = Id)       
           )

#### latent variable model ####
library(lava)
e.lvm <- estimate(lvm(Y ~ G), data = data.df)
extractData(e.lvm)
extractData(e.lvm, design.matrix = TRUE)

Description

Find all new links between variables (adapted from lava::modelsearch).

Usage

findNewLink(object, ...)

## S3 method for class 'lvm'
findNewLink(
  object,
  data = NULL,
  type = "both",
  exclude.var = NULL,
  rm.latent_latent = FALSE,
  rm.endo_endo = FALSE,
  rm.latent_endo = FALSE,
  output = "names",
  ...
)

Arguments

Value

A list containing:

• M.links: a matrix with two columns indicating (by name or position) the exogenous and endogenous variable corresponding to each link.

• links: the name of the additional possible links

• directional: a logical vector indicating for each link whether the link is unidirectional (TRUE, i.e. regression link) or bidirectional (FALSE, i.e. covariance link).

Examples

library(lava)

m <- lvm()
regression(m) <- c(y1,y2,y3)~u
categorical(m,labels=c("M","F","MF")) <- ~X1
findNewLink(m, rm.endo = FALSE)
findNewLink(m, rm.endo = TRUE)
findNewLink(m, exclude.var = "X1")

regression(m) <- u~x1+x2
latent(m) <- ~u

findNewLink(m, rm.endo = FALSE)
findNewLink(m, rm.endo = TRUE)
findNewLink(m, rm.endo = TRUE, output = "index")
findNewLink(m, type = "covariance")
findNewLink(m, type = "regression")

Description

Estimate LVM with weights.

Usage

gaussian_weight.estimate.hook(x, data, estimator, ...)

gaussian_weight_method.lvm

gaussian_weight_logLik.lvm(object, type = "cond", p, data, weights, ...)

gaussian_weight_objective.lvm(x, ...)

gaussian_weight_score.lvm(
  x,
  data,
  p,
  S,
  n,
  mu = NULL,
  weights = NULL,
  debug = FALSE,
  reindex = FALSE,
  mean = TRUE,
  constrain = TRUE,
  indiv = FALSE,
  ...
)

gaussian_weight_gradient.lvm(...)

gaussian_weight_hessian.lvm(x, p, n, weights = NULL, ...)

Arguments

Format

An object of class character of length 1.

Examples

#### linear regression with weights ####

## data
df <- data.frame(Y = c(1,2,2,1,2),
                 X = c(1,1,2,2,2),
                 missing = c(0,0,0,0,1),
                 weights = c(1,1,2,1,NA))

## using lm
e.lm.GS <- lm(Y~X, data = df)
e.lm.test <- lm(Y~X, data = df[df$missing==0,], weights = df[df$missing==0,"weights"])

## using lvm
m <- lvm(Y~X)
e.GS <- estimate(m, df)
## e.lava.test <- estimate(m, df[df$missing==0,], weights = df[df$missing==0,"weights"])
## warnings!!
e.test <- estimate(m, data = df[df$missing==0,],
                   weights = df[df$missing==0,"weights"],
                   estimator = "gaussian_weight")

Description

Extract the links that have been found relevant by modelsearch2.

Usage

getNewLink(object, step)

## S3 method for class 'modelsearch2'
getNewLink(object, step = 1:nStep(object))

Arguments

Value

A character vector.

Examples

## Not run: 
mSim <- lvm(Y~G+X1+X2)
addvar(mSim) <- ~Z1+Z2+Z3+Z4+Z5+Z6

set.seed(10)
df.data <- lava::sim(mSim, 1e2)

mBase <- lvm(Y~G)
addvar(mBase) <- ~X1+X2+Z1+Z2+Z3+Z4+Z5+Z6
e.lvm <- estimate(mBase, data = df.data)
res <- modelsearch2(e.lvm, method.p.adjust = "holm")
getNewLink(res)

## End(Not run)

Description

Extract the model that has been retained by modelsearch2.

Usage

getNewModel(object, step)

## S3 method for class 'modelsearch2'
getNewModel(object, step = nStep(object))

Arguments

Value

A lvmfit object.

Examples

## Not run: 
mSim <- lvm(Y~G+X1+X2)
addvar(mSim) <- ~Z1+Z2+Z3+Z4+Z5+Z6

set.seed(10)
df.data <- lava::sim(mSim, 1e2)

mBase <- lvm(Y~G)
addvar(mBase) <- ~X1+X2+Z1+Z2+Z3+Z4+Z5+Z6
e.lvm <- estimate(mBase, data = df.data)
res <- modelsearch2(e.lvm, method.p.adjust = "holm")
getNewModel(res)

## End(Not run)

Description

Extract one step from the sequential procedure.

Usage

getStep(object, step, slot)

## S3 method for class 'modelsearch2'
getStep(object, step = nStep(object), slot = NULL)

Arguments

Examples

## Not run: 
mSim <- lvm(Y~G+X1+X2)
addvar(mSim) <- ~Z1+Z2+Z3+Z4+Z5+Z6
df.data <- lava::sim(mSim, 1e2)

mBase <- lvm(Y~G)
addvar(mBase) <- ~X1+X2+Z1+Z2+Z3+Z4+Z5+Z6
e.lvm <- estimate(mBase, data = df.data)
res <- modelsearch2(e.lvm, method.p.adjust = "holm")

getStep(res)
getStep(res, slot = "sequenceTest")
getStep(res, step = 1)

## End(Not run)

Description

Reconstruct the residual variance-covariance matrix from a latent variable model. It is similar to nlme::getVarCov but with small sample correction.

Usage

getVarCov2(object, ...)

## S3 method for class 'lvmfit'
getVarCov2(object, ssc = lava.options()$ssc, ...)

## S3 method for class 'lvmfit2'
getVarCov2(object, ...)

Arguments

Details

When argument object is a lvmfit object, the method first calls estimate2 and then extract the residuals.

Value

A matrix with as many rows and column as the number of endogenous variables

Examples

#### simulate data ####
set.seed(10)
n <- 101

Y1 <- rnorm(n, mean = 0)
Y2 <- rnorm(n, mean = 0.3)
Id <- findInterval(runif(n), seq(0.1,1,0.1))
data.df <- rbind(data.frame(Y=Y1,G="1",Id = Id),
           data.frame(Y=Y2,G="2",Id = Id)
           )

#### latent variable models ####
library(lava)
e.lvm <- estimate(lvm(Y ~ G), data = data.df)
getVarCov2(e.lvm)

Description

Test linear hypotheses on coefficients from a latent variable models with small sample corrections.

Usage

glht2(object, ...)

## S3 method for class 'lvmfit'
glht2(
  object,
  linfct,
  rhs = NULL,
  robust = FALSE,
  cluster = NULL,
  ssc = lava.options()$ssc,
  df = lava.options()$df,
  ...
)

## S3 method for class 'lvmfit2'
glht2(object, linfct, rhs = NULL, robust = FALSE, cluster = NULL, ...)

## S3 method for class 'mmm'
glht2(object, linfct, rhs = 0, robust = FALSE, cluster = NULL, ...)

## S3 method for class 'lvmfit2'
glht(model, linfct, rhs = NULL, robust = FALSE, cluster = NULL, ...)

Arguments

Details

Whenever the argument linfct is not a matrix, it is passed to the function createContrast to generate the contrast matrix and, if not specified, rhs.

Since only one degree of freedom can be specify in a glht object and it must be an integer, the degree of freedom of the denominator of an F test simultaneously testing all hypotheses is retained, after rounding.

Argument rhs and null are equivalent. This redondance enable compatibility between lava::compare, compare2, multcomp::glht, and glht2.

Value

A glht object.

See Also

createContrast to create contrast matrices.
estimate2 to pre-compute quantities for the small sample correction.

Examples

library(multcomp)

## Simulate data
mSim <- lvm(c(Y1,Y2,Y3)~ beta * eta, Z1 ~ E, Z2 ~ E, Age[40:5]~1)
latent(mSim) <- "eta"
set.seed(10)
n <- 1e2

df.data <- lava::sim(mSim, n, latent = FALSE, p = c(beta = 1))

#### Inference on a single model ####
e.lvm <- estimate(lvm(Y1~E), data = df.data)
summary(glht2(e.lvm, linfct = c("Y1~E + Y1","Y1")))

#### Inference on separate models ####
## fit separate models
lvmX <- estimate(lvm(Z1 ~ E), data = df.data)
lvmY <- estimate(lvm(Z2 ~ E + Age), data = df.data)
lvmZ <- estimate(lvm(c(Y1,Y2,Y3) ~ eta, eta ~ E), 
                 data = df.data)

#### create mmm object #### 
e.mmm <- mmm(X = lvmX, Y = lvmY, Z = lvmZ)

#### create contrast matrix ####
resC <- createContrast(e.mmm, linfct = "E")

#### adjust for multiple comparisons ####
e.glht2 <- glht2(e.mmm, linfct = c(X="E"), df = FALSE)
summary(e.glht2)

Description

Extract the hessian from a latent variable model, with small sample correction

Usage

hessian2(object, indiv, cluster, as.lava, ...)

## S3 method for class 'lvmfit'
hessian2(
  object,
  indiv = FALSE,
  cluster = NULL,
  as.lava = TRUE,
  ssc = lava.options()$ssc,
  ...
)

## S3 method for class 'lvmfit2'
hessian2(object, indiv = FALSE, cluster = NULL, as.lava = TRUE, ...)

Arguments

Details

When argument object is a lvmfit object, the method first calls estimate2 and then extract the hessian.

Value

An array containing the second derivative of the likelihood relative to each sample (dim 3) and each pair of model coefficients (dim 1,2).

See Also

estimate2 to obtain lvmfit2 objects.

Examples

#### simulate data ####
n <- 5e1
p <- 3
X.name <- paste0("X",1:p)
link.lvm <- paste0("Y~",X.name)
formula.lvm <- as.formula(paste0("Y~",paste0(X.name,collapse="+")))

m <- lvm(formula.lvm)
distribution(m,~Id) <- Sequence.lvm(0)
set.seed(10)
d <- lava::sim(m,n)

#### latent variable models ####
e.lvm <- estimate(lvm(formula.lvm),data=d)
hessian2(e.lvm)

Description

Extract the influence function from a latent variable model. It is similar to lava::iid but with small sample correction.

Usage

iid2(object, ...)

## S3 method for class 'lvmfit'
iid2(
  object,
  robust = TRUE,
  cluster = NULL,
  as.lava = TRUE,
  ssc = lava.options()$ssc,
  ...
)

## S3 method for class 'lvmfit2'
iid2(object, robust = TRUE, cluster = NULL, as.lava = TRUE, ...)

## S3 method for class 'lvmfit2'
iid(x, robust = TRUE, cluster = NULL, as.lava = TRUE, ...)

Arguments

Details

When argument object is a lvmfit object, the method first calls estimate2 and then extract the variance-covariance matrix.

Value

A matrix containing the 1st order influence function relative to each sample (in rows) and each model coefficient (in columns).

See Also

estimate2 to obtain lvmfit2 objects.

Examples

#### simulate data ####
n <- 5e1
p <- 3
X.name <- paste0("X",1:p)
link.lvm <- paste0("Y~",X.name)
formula.lvm <- as.formula(paste0("Y~",paste0(X.name,collapse="+")))

m <- lvm(formula.lvm)
distribution(m,~Id) <- Sequence.lvm(0)
set.seed(10)
d <- sim(m,n)

#### latent variable model ####
e.lvm <- estimate(lvm(formula.lvm),data=d)
iid.tempo <- iid2(e.lvm)

Description

Extract the i.i.d. decomposition and display it along with the corresponding coefficient.

Usage

iid2plot(object, param)

Arguments

Description

Extract iid decomposition (i.e. influence function) from model object.

Usage

iidJack(object, ...)

## Default S3 method:
iidJack(
  object,
  data = NULL,
  grouping = NULL,
  cpus = 1,
  keep.warnings = TRUE,
  keep.error = TRUE,
  cl = NULL,
  trace = TRUE,
  ...
)

Arguments

Value

A matrix with in row the samples and in columns the parameters.

Examples

n <- 20
set.seed(10)
mSim <- lvm(c(Y1,Y2,Y3,Y4,Y5) ~ 1*eta)
latent(mSim) <- ~eta
categorical(mSim, K=2) <- ~G
transform(mSim, Id ~ eta) <- function(x){1:NROW(x)}
dW <- lava::sim(mSim, n, latent = FALSE)

#### LVM ####
## Not run: 
m1 <- lvm(c(Y1,Y2,Y3,Y4,Y5) ~ 1*eta)
latent(m1) <- ~eta
regression(m1) <- eta ~ G
e <- estimate(m1, data = dW)

iid1 <- iidJack(e)
iid2 <- iid(e)
attr(iid2, "bread") <- NULL

apply(iid1,2,sd)
apply(iid2,2,sd)
quantile(iid2 - iid1)

## End(Not run)

Description

Extract the expected information matrix from a latent variable model. Similar to lava::information but with small sample correction.

Usage

information2(object, as.lava, ssc, ...)

## S3 method for class 'lvmfit'
information2(object, as.lava = TRUE, ssc = lava.options()$ssc, ...)

## S3 method for class 'lvmfit2'
information2(object, as.lava = TRUE, ...)

## S3 method for class 'lvmfit2'
information(x, ...)

Arguments

Details

When argument object is a lvmfit object, the method first calls estimate2 and then extract the information matrix.

Value

A matrix with as many rows and columns as the number of coefficients.

See Also

estimate2 to obtain lvmfit2 objects.

Examples

#### simulate data ####
n <- 5e1
p <- 3
X.name <- paste0("X",1:p)
link.lvm <- paste0("Y~",X.name)
formula.lvm <- as.formula(paste0("Y~",paste0(X.name,collapse="+")))

m <- lvm(formula.lvm)
distribution(m,~Id) <- Sequence.lvm(0)
set.seed(10)
d <- lava::sim(m,n)

#### linear models ####
e.lm <- lm(formula.lvm,data=d)

#### latent variable models ####
e.lvm <- estimate(lvm(formula.lvm),data=d)
information(e.lvm)
information2(e.lvm)
information2(e.lvm)[1:4,1:4] -  solve(vcov(e.lm))

Description

Convert var1 and var2 from formula or covariance to character.

Usage

initVarLink(
  var1,
  var2,
  rep.var1 = FALSE,
  format = "list",
  Slink = c(lava.options()$symbols[1], "~"),
  Scov = lava.options()$symbols[2]
)

initVarLinks(var1, format = "list", ...)

Arguments

Value

See argument format.

Examples

initVarLink(y ~ x1)
initVarLink("y ~ x1")
initVarLink(y ~ x1 + x2)
initVarLink("y ~ x1 + x2")
initVarLink(y ~ x1 + x2, rep.var1 = TRUE)
initVarLink(y ~ x1 + x2, rep.var1 = TRUE, format = "formula")
initVarLink(y ~ x1 + x2, rep.var1 = TRUE, format = "txt.formula")
initVarLink("y", "x1", format = "formula")

initVarLink("y ~ x1:0|1")

initVarLinks(y ~ x1)
initVarLinks("y ~ x1")
initVarLinks(c("y ~ x1","y~ x2"))
initVarLinks(c(y ~ x1,y ~ x2))
initVarLinks(c("y ~ x1","y ~ x2"), format = "formula")
initVarLinks(c(y ~ x1,y ~ x2), format = "formula")
initVarLinks(c("y ~ x1","y~ x2"), format = "txt.formula")
initVarLinks(c(y ~ x1,y ~ x2), format = "txt.formula")

Description

Consider a univariate random variable X, two multivariate random variables Y and Z, and t1 and t2 two real numbers. This function can compute either P[|X|>t1,|X]>|Y1|,...,|X]>|Yp|] if zmin is not specified, P[|Z1|<t2,...,|Zq|<t2,|X|>t1,|X]>|Y1|,...,|X]>|Yp|] if zmin is specified.

Usage

intDensTri(
  mu,
  Sigma,
  df,
  n,
  x.min,
  z.max = NULL,
  type = "double",
  proba.min = 1e-06,
  prune = NULL,
  distribution = "pmvnorm"
)

Arguments

Details

Argument type:

• \"raw\": mesh with points inside the domain

• \"double\": mesh with points outside the domain

• \"fine\": mesh with points inside the domain plus additional rectangles trying to fill the missing domain.

Argument Sigma and mu: define the mean and variance-covariance of the random variables X, Y, Z (in this order). The length of the argument z.max is used to define the dimension of Z. The dimension of X is always 1.

Value

A numeric.

Examples

library(mvtnorm)

p <- 2
Sigma <- diag(p)
mu <- rep(0, p)

## bivariate normal distribution
z2 <- qmvt(0.975, mean = mu, sigma = Sigma, df = 1e3)$quantile

# compute integral
intDensTri(mu = mu, Sigma = Sigma, n=5, x.min=0, type = "fine")$value-1/2
intDensTri(mu = mu, Sigma = Sigma, n=30, x.min=0, type = "raw")$value-1/2
intDensTri(mu = mu, Sigma = Sigma, n=50, x.min=0, type = "raw")$value-1/2

intDensTri(mu = mu, Sigma = Sigma, df = 5, n=5, x.min=0, distribution = "pmvt")$value-1/2
res <- intDensTri(mu = mu, Sigma = Sigma, df = 5, n=10, x.min=0, distribution = "pmvt")
res$value-1/2
ggplot2::autoplot(res)

## trivariate normal distribution
## Not run: 
p <- 3
Sigma <- diag(p)
mu <- rep(0, p)

res2 <- intDensTri(mu = mu, Sigma = Sigma, n=5, x.min = 0, z.max = 10)
ggplot2::autoplot(res2)
ggplot2::autoplot(res2, coord.plot = c("x","z1"))
res2

## End(Not run)

#### when the distribution is far from 0
## Not run: 
eq1 <- intDensTri(mu = c(10,0), Sigma = diag(1,2), 
                  x.min = 2, n=10)
eq1$value-1
ggplot2::autoplot(eq1)

eq2 <- intDensTri(mu = c(10,0,0), Sigma = diag(1,3),
                  x.min=2, z.max = 10, type = "raw",
                  n=10)
ggplot2::autoplot(eq2, coord.plot = c("y1","z1"))
eq2$value-1

## more variables
p <- 5
Sigma <- diag(p)
mu <- rep(0, p)

res2 <- intDensTri(mu = mu, Sigma = Sigma, n=5, x.min = 1, z.max = c(2,2))
res2$grid

## End(Not run)

Description

The package contains three main functionalities:

• compare2, summary2: Wald tests/robust Wald tests/F-tests/robust F-tests with improved control of the type 1 error in small samples.

• glht2: adjustment for multiple comparisons when doing inference for multiple latent variable models.

• modelsearch2: searching for local dependencies with adjustment for multiple comparisons.

It contains other useful functions such as:

• calibrateType1: simulation study of the type 1 error of Wald tests.

• createContrast: user-friendly function generating a contrast matrix.

• getVarCov2: reconstruct the conditional variance covariance matrix.

• iidJack: extract the jackknife iid decomposition.

Details

The latent variable models (LVM) considered in this package can be written
as a measurement model:

$Y_i = \nu + \eta_i \Lambda + X_i K + \epsilon_i$

and a structural model:

$\eta_i = \alpha + \eta_i B + X_i \Gamma + \zeta_i$

where $\Sigma$ is the variance covariance matrix of the residuals $\epsilon$,
and $\Psi$ is the variance covariance matrix of the residuals $\zeta$.

The corresponding conditional mean is:

$\mu_i(\theta) = E[Y_i|X_i] = \nu + (\alpha + X_i \Gamma) (1-B)^{-1} \Lambda + X_i K$

$\Omega(\theta) = Var[Y_i|X_i] = \Lambda^{t} (1-B)^{-t} \Psi (1-B)^{-1} \Lambda + \Sigma$

The package aims to provides tool for testing linear hypotheses on the model coefficients $\nu$, $\Lambda$, $K$, $\Sigma$, $\alpha$, $B$, $\Gamma$, $\Psi$. Searching for local dependency enable to test whether the proposed model is too simplistic and if so to identify which additional coefficients should be added to the model.

Limitations

'lavaSearch2' has been design for Gaussian latent variable models. This means that it may not work / give valid results:

• in presence of censored or binary outcomes.

• with stratified models (i.e. object of class multigroup).

Description

Extract leverage values from a latent variable model, with small sample correction.

Usage

leverage2(object, format, ssc, ...)

## S3 method for class 'lvmfit'
leverage2(object, format = "wide", ssc = lava.options()$ssc, ...)

## S3 method for class 'lvmfit2'
leverage2(object, format = "wide", ...)

Arguments

Details

The leverage are defined as the partial derivative of the fitted values with respect to the observations.

$leverage_i = \frac{\partial \hat{Y}_i}{\partial Y_i}$

See Wei et al. (1998).

When argument object is a lvmfit object, the method first calls estimate2 and then extract the leverage.

Value

a matrix containing the leverage relative to each sample (in rows) and each endogenous variable (in column).

References

Bo-Cheng Wei et al., Generalized Leverage and its applications (1998), Scandinavian Journal of Statistics 25:1:25-37.

See Also

estimate2 to obtain lvmfit2 objects.

Examples

#### simulate data ####
set.seed(10)
m <- lvm(Y1~eta,Y2~eta,Y3~eta)
latent(m) <- ~eta
d <- lava::sim(m,20, latent = FALSE)

#### latent variable models ####
e.lvm <- estimate(m, data = d)
leverage2(e.lvm)

Description

Compute the power of a matrix.

Usage

matrixPower(object, power, symmetric, tol = 1e-12, print.warning = TRUE)

Arguments

Value

A matrix.

Examples

## symmetric matrix
set.seed(10)
M <- matrix(rnorm(20*6),20,6)
Sigma <- var(M)
Sigma.half <- matrixPower(Sigma, power = 1/2, symmetric = TRUE)
round(Sigma.half %*% Sigma.half - Sigma,5)

iSigma <- matrixPower(Sigma, power = -1, symmetric = TRUE)
round(iSigma %*% Sigma,5)

iSigma.half <- matrixPower(Sigma, power = -1/2, symmetric = TRUE)
round(iSigma.half %*% iSigma.half - iSigma,5)

## non symmetric matrix
set.seed(10)
M <- matrix(abs(rnorm(9)), 3, 3) + diag(1,3,3)
M-t(M)

iM <- matrixPower(M, power = -1, symmetric = FALSE)
round(iM %*% M,5)

iM.half <- matrixPower(M, power = -1/2, symmetric = FALSE)
round(iM.half %*% iM.half %*% M,5)

Description

Procedure adding relationship between variables that are supported by the data.

Usage

modelsearch2(
  object,
  link,
  data,
  method.p.adjust,
  method.maxdist,
  n.sample,
  na.omit,
  alpha,
  nStep,
  trace,
  cpus
)

## S3 method for class 'lvmfit'
modelsearch2(
  object,
  link = NULL,
  data = NULL,
  method.p.adjust = "fastmax",
  method.maxdist = "approximate",
  n.sample = 1e+05,
  na.omit = TRUE,
  alpha = 0.05,
  nStep = NULL,
  trace = TRUE,
  cpus = 1
)

Arguments

Details

method.p.adjust = "max" computes the p-values based on the distribution of the max statistic. This max statistic is the max of the square root of the score statistic. The p-value are computed integrating the multivariate normal distribution.

method.p.adjust = "fastmax" only compute the p-value for the largest statistic. It is faster than "max" and lead to identical results.

method.p.adjust = "gof" keep adding links until the chi-squared test (of correct specification of the covariance matrix) is no longer significant.

Value

A list containing:

• sequenceTest: the sequence of test that has been performed.

• sequenceModel: the sequence of models that has been obtained.

• sequenceQuantile: the sequence of rejection threshold. Optional.

• sequenceIID: the influence functions relative to each test. Optional.

• sequenceSigma: the covariance matrix relative to each test. Optional.

• initialModel: the model before the sequential search.

• statistic: the argument statistic.

• method.p.adjust: the argument method.p.adjust.

• alpha: [numeric 0-1] the significance cutoff for the p-values.

• cv: whether the procedure has converged.

Examples

## simulate data 
mSim <- lvm()
regression(mSim) <- c(y1,y2,y3,y4)~u
regression(mSim) <- u~x1+x2
categorical(mSim,labels=c("A","B","C")) <- "x2"
latent(mSim) <- ~u
covariance(mSim) <- y1~y2
transform(mSim, Id~u) <- function(x){1:NROW(x)}

set.seed(10)
df.data <- lava::sim(mSim, n = 1e2, latent = FALSE)

## only identifiable extensions
m <- lvm(c(y1,y2,y3,y4)~u)
latent(m) <- ~u
addvar(m) <- ~x1+x2

e <- estimate(m, df.data)

## Not run: 
resSearch <- modelsearch(e)
resSearch

resSearch2 <- modelsearch2(e, nStep = 2)
resSearch2

## End(Not run)


## some extensions are not identifiable
m <- lvm(c(y1,y2,y3)~u)
latent(m) <- ~u
addvar(m) <- ~x1+x2 

e <- estimate(m, df.data)

## Not run: 
resSearch <- modelsearch(e)
resSearch
resSearch2 <- modelsearch2(e)
resSearch2

## End(Not run)

## for instance
mNI <- lvm(c(y1,y2,y3)~u)
latent(mNI) <- ~u
covariance(mNI) <- y1~y2
## estimate(mNI, data = df.data)
## does not converge

Description

Extract the effective sample size, i.e. sample size minus the loss in degrees of freedom caused by the estimation of the parameters.

Usage

nobs2(object, ssc, ...)

## S3 method for class 'lvmfit'
nobs2(object, ssc = lava.options()$ssc, ...)

## S3 method for class 'lvmfit2'
nobs2(object, ...)