## ----setup,include=FALSE------------------------------------------------------
library("lava")
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
  )
has_mets <- lava:::versioncheck('mets', 1)
has_geepack <- lava:::versioncheck('geepack', 1)

## ----sim_model----------------------------------------------------------------
m <- lvm() |>
  regression(y1 ~ x1 + a + w) |>
  regression(y2 ~ x2 + a + w) |>
  regression(y3 ~ x3 + a + w) |>
  regression(y4 ~ x4 + a + w) |>
  regression(a ~ w) |>
  distribution(~ y1 + y2 + y3 + y4 + a, value = binomial.lvm()) |>
  distribution(~id, value = Sequence.lvm(integer = TRUE))

## ----simulate-----------------------------------------------------------------
n <- 4e2
dw <- sim(m, n, seed = 1) |>
  transform(y3 = y3 * ifelse(id > n / 2, NA, 1))
Print(dw)
## Data in long format
dl <- reshape(dw,
        varying = list(paste0("y",1:4),
                       paste0("x",1:4)),
        v.names=c("y", "x"), direction="long") |>
  na.omit()
dl <- dl[order(dl$id), ]
## dl <- mets::fast.reshape(dw, varying = c("y", "x")) |> na.omit()
Print(dl)

## ----estimate.syntax, eval=FALSE----------------------------------------------
# estimate(x=, ...)
# estimate(coef=, IF=, ...)
# estimate(coef=, vcov=, ...)

## ----inp1---------------------------------------------------------------------
inp <- as.matrix(dw[, c("y1", "y2")])
e <- estimate(inp[, 1, drop = FALSE], type="mean") 
class(e)
e

## ----ic1----------------------------------------------------------------------
IC(e) |> Print()

## ----inp2---------------------------------------------------------------------
estimate(inp)

## ----mlm----------------------------------------------------------------------
e <- lm(cbind(y1, y2) ~ 1, data = dw) |>
  estimate()
IC(e) |> head()

## ----estimatemethods----------------------------------------------------------
summary(e)
## extract parameter coefficients
coef(e)
## ## Asymptotic (robust) variance estimate
vcov(e)
## Matrix with estimates and confidence limits
estimate(e, level = 0.99) |> parameter()
## Influence curve
IC(e) |> head()
## Join estimates
e + e # Same as merge(e,e)

## ----glm----------------------------------------------------------------------
g <- glm(y1 ~ a + x1, data = dw, family = binomial)
estimate(g)

## ----glm.std------------------------------------------------------------------
estimate(g, robust = FALSE)

## ----ifglm--------------------------------------------------------------------
IC(g) |> head()

## ----ordreg-------------------------------------------------------------------
ordreg(y1 ~ a + x1, dw, family=binomial(logit)) |> estimate()

## ----mets---------------------------------------------------------------------
library("survival")
data(pbc, package="survival")

## ----phreg, eval=has_mets-----------------------------------------------------
fit.phreg <- mets::phreg(Surv(time, status > 0) ~ age + sex, data = pbc)
fit.phreg
IC(fit.phreg) |> head()

## ----phreg-baseline, eval=has_mets--------------------------------------------
baseline <- function(object, time, ...) {
  ic <- mets::IC(object, baseline = TRUE, time = time, ...)
  est <- mets::predictCumhaz(object$cumhaz, new.time = time)[1, 2]
  estimate(NULL, coef = est, IC = ic, labels = paste0("chaz:", time))
}
tt <- 2000
baseline(fit.phreg, tt)

## ----survreg, eval=has_mets---------------------------------------------------
survival::survreg(Surv(time, status > 0) ~ age + sex, data = pbc, dist="weibull") |>
  estimate()

## ----semfit, eval=has_mets----------------------------------------------------
sem <- lvm(y1 + y2 ~ 1 * u + w) |>
  latent(~ u) |>
  ordinal(K=2, ~ y1 + y2)
semfit <- estimate(sem, data = dw)

## Robust standard errors
estimate(semfit)

## ----quantiles----------------------------------------------------------------
eq <- estimate(dw[, c("w", "x1")], type = "quantile", probs = c(0.25, 0.5, 0.75))
eq
IC(eq) |> head()

## ----glmmarg------------------------------------------------------------------
g1 <- glm(y1 ~ a, family=binomial, data=dw)
g2 <- glm(y2 ~ a, family=binomial, data=dw)
e <- merge(g1, g2)
summary(e)

## ----hypo1--------------------------------------------------------------------
estimate(e, cbind(0,1,0,-1), null=0)

## ----glmmargmis---------------------------------------------------------------
g2 <- glm(y2 ~ 1, family = binomial, data = dw)
summary(g2)
dwc <- na.omit(dw) 
g3 <- glm(y3 ~ 1, family = binomial, data = dwc)
summary(g3)

e2 <- estimate(g2, id = dw$id)
e3 <- estimate(g3, id = "id", data=dwc)

merge(e2,e3) |> IC() |> Print()
vcov(e2 + e3)
## Same marginals as
list(vcov(e2), vcov(e3))

## ----merge--------------------------------------------------------------------
merge(e2, e3, id = list(dw$id, dwc$id))

## ----estimatenoid-------------------------------------------------------------
estimate(g2) |>
  IC() |> head()
vcov(estimate(g2) + estimate(g3))
(estimate(g2) + estimate(g3)) |>
  (rownames %++% head %++% IC)()

## ----merge_idnull-------------------------------------------------------------
merge(g1, g2, id = NULL) |> (Print %++% IC)()
merge(g1, g2, id = NULL) |> vcov()

## ----mergekeep----------------------------------------------------------------
merge(g1, g2, keep = c("(Intercept)", "(Intercept).1"))

## ----merge2-------------------------------------------------------------------
merge(g1,g2, keep=c(1, 3))

## ----merge3-------------------------------------------------------------------
merge(g1, g2, keep = "cept", regex = TRUE)
merge(g1, g2, keep = c("\\)$", "^a$"), regex = TRUE, ignore.case = TRUE)

## ----merge4-------------------------------------------------------------------
merge(g1, g2, labels = c("a", "b", "c")) |> estimate(keep = c("a", "c"))
merge(g1, g2,
      labels = c("a", "b", "c"),
      keep = c("a", "c")
)
estimate(g1, labels=c("a", "b"))

## ----merge5-------------------------------------------------------------------
merge(g1, g2, subset="(Intercept)")

## ----cluster1-----------------------------------------------------------------
g0 <- glm(y ~ a + w + x, data = dl, family = binomial())

## ----cluster2-----------------------------------------------------------------
estimate(g0, id=dl$id)

## ----geepack, eval=has_geepack------------------------------------------------
gee0 <- geepack::geeglm(y ~ a + w + x, data = dl, id = dl$id, family=binomial)
summary(gee0)

## ----aggregate----------------------------------------------------------------
set.seed(1)
y <- cbind(rnorm(1e5))
N <- 2e2 ## Number of aggregated groups, the number of observations in the new IF
id <- foldr(nrow(y), N, list=FALSE)
Print(cbind(table(id)))

## Aggregated IF
e <- estimate(cbind(y), id = id) 
object.size(e)
e

## ----delta1-------------------------------------------------------------------
estimate(g1, sum)
estimate(g1, function(p) list(a = sum(p))) # named list
## Multiple parameters
estimate(g1, function(x) c(x, x[1] + exp(x[2]), inv = 1 / x[2]))
estimate(g1, exp)         

## ----cov----------------------------------------------------------------------
Cov <- function(x, y, ...) {
  est <- mean(x * y)-mean(x)*mean(y)
    estimate(
      coef = est,
      IC = (x - mean(x)) * (y - mean(y)) - est,
      ...
    )
}
with(dw, Cov(x1, x2))

## ----cov2---------------------------------------------------------------------
e1 <- lm(cbind(x1, x2, x1 * x2) ~ 1, data = dw) |>
  estimate(labels = c("Ex1", "Ex2", "Ex1x2"))
e1
estimate(e1, function(x) c(x, cov=with(as.list(x), Ex1x2 - Ex2* Ex1)))

## ----rho----------------------------------------------------------------------
e2 <- with(dw, Cov(x1, x2, labels = "c", id = id) +
               Cov(x1, x1, labels = "v1", id = id) +
               Cov(x2, x2, labels = "v2", id = id))
rho <- estimate(e2, function(x) list(rho = x[1] / (x[2] * x[3])^.5))
rho

## ----tanh---------------------------------------------------------------------
estimate(rho, atanh, back.transform = tanh)

## -----------------------------------------------------------------------------
g <- lapply(
  list(y1 ~ a, y2 ~ a, y3 ~ a), #, y4 ~ a+x4),
  function(f) glm(f, family = binomial, data = dw)
)
gg <- Reduce(merge, g)
gg

## ----contrast1----------------------------------------------------------------
B <- cbind(0,1, 0,-1, 0,0)
estimate(gg, B)

## ----estcontrast1-------------------------------------------------------------
estimate(gg, B, null=1)

## ----estcontrast2-------------------------------------------------------------
B <- rbind(cbind(0,1, 0,-1, 0,0),
           cbind(0,1, 0,0, 0,-1))
estimate(gg, B)

## -----------------------------------------------------------------------------
estimate(gg, a + a.1, 2*a - a.2, a, null=c(2,1,1))

## ----contr--------------------------------------------------------------------
contr(list(1, c(1, 2), c(1, 4)), n = 5)

## ----pairwise.diff------------------------------------------------------------
pairwise.diff(3)
estimate(gg, pairwise.diff(3), null=c(1,1,1), use=c(2,4,6))

## ----pcorrect, eval=has_mets--------------------------------------------------
gg0 <- estimate(gg, use="^a", regex=TRUE, null=rep(.8, 3))
p.correct(gg0)

## ----closedtesting, eval=has_mets---------------------------------------------
closed.testing(gg0)

## ----estpred------------------------------------------------------------------
g <- glm(y1 ~ a + x1 + w, data=dw, family=binomial)
pr <- function(p, data, ...)
  with(data, expit(p[1] + p["a"] + p["x1"]*x1 + p["w"]*w))
pr(coef(g), dw) |> head()

## ----average------------------------------------------------------------------
id <- foldr(NROW(dw), 100, list=FALSE)
ea <- estimate(g, pr, average=TRUE, id=id)
ea
IC(ea) |> head()

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