## ----include=FALSE------------------------------------------------------------
mets <- lava:::versioncheck('mets', 1)


## ----plotload, results="hide",message=FALSE,warning=FALSE---------------------
library("lava")
m0 <- lvm(y1+y2 ~ x, y1 ~~ y2)
edgelabels(m0, y1 + y2 ~ x) <- c(expression(b[1]), expression(b[2]))
edgelabels(m0, y1 ~ y2) <- "r"
plot(m0, layoutType="circo")


## ----m0-----------------------------------------------------------------------
m0 <- lvm() |>
  covariance(y1 ~ y2, value='r') |>
  regression(y1 + y2 ~ x)


## ----coef---------------------------------------------------------------------
coef(m0, labels=TRUE)


## ----sim----------------------------------------------------------------------
d <- sim(m0, 500, p=c(r=0.9), seed=1)
head(d)


## ----defcens------------------------------------------------------------------
cens1 <- function(threshold,type='right') {
  function(x) {
    x <- unlist(x)
    if (type=='left')
      return( survival::Surv(pmax(x,threshold), x>=threshold, type='left') )
      return ( survival::Surv(pmin(x,threshold), x<=threshold) )
  }
}

m0 <-
  transform(m0, s1 ~ y1, cens1(-2, 'left')) |>
  transform(s2 ~ y2, cens1(2,  'right'))


## ----sim2---------------------------------------------------------------------
d <- sim(m0, 500, p=c(r=0.9), seed=1)
head(d)


## ----est1---------------------------------------------------------------------
m <- lvm() |>
     regression(y1 + y2 ~ x) |>
     covariance(y1 ~ y2)

e <- estimate(m, data=d)
e


## ----delta--------------------------------------------------------------------
estimate(e, function(p) p['y1~~y2']/(p['y1~~y1']*p['y2~~y2'])^.5)


## ----correlation--------------------------------------------------------------
correlation(e)


## -----------------------------------------------------------------------------
estimate(e, function(p) atanh(p['y1~~y2']/(p['y1~~y1']*p['y2~~y2'])^.5)) |>
  summary(transform=tanh)


## ----constraints--------------------------------------------------------------
m2 <- m |>
    parameter(~ l1 + l2 + z) |>
    variance(~ y1 + y2, value=c('v1','v2')) |>
    covariance(y1 ~ y2, value='c') |>
    constrain(v1 ~ l1, fun=exp) |>
    constrain(v2 ~ l2, fun=exp) |>
    constrain(c ~ z+l1+l2, fun=function(x) tanh(x[1])*sqrt(exp(x[2])*exp(x[3])))


## ----estconstraints-----------------------------------------------------------
e2 <- estimate(m2, d)
e2


## ----deltaconstraints---------------------------------------------------------
estimate(e2, "z") |> labels("rho") |> summary(transform=tanh)


## ----constraints2-------------------------------------------------------------
m2 <- lvm() |>
  regression(y1 + y2 ~ x) |>
  covariance(y1 ~ y2, constrain=TRUE, rname='z')

e2 <- estimate(m2, data=d)
e2


## ----e2backtransform----------------------------------------------------------
estimate(e2, 'z') |> summary(transform=tanh)


## ----profileci, cache=TRUE----------------------------------------------------
tanh(confint(e2, 'z', profile=TRUE))


## ----bootstrap, cache=TRUE----------------------------------------------------
set.seed(1)
b <- bootstrap(e2, data=d, R=50, mc.cores=1)
b


## ----cache=TRUE---------------------------------------------------------------
quantile(tanh(b$coef[,'z']), c(.025,.975))
















## -----------------------------------------------------------------------------
sessionInfo()