misaem (version 0.9.0)

miss.saem: miss.saem

Description

This function uses algorithm SAEM to fit the logistic regression model with missing data.

Usage

miss.saem(X.obs, y, pos_var = 1:ncol(X.obs), maxruns = 500,
  tol_em = 1e-07, nmcmc = 2, tau = 1, k1 = 50, seed = 200,
  print_iter = TRUE, var_cal = FALSE, ll_obs_cal = FALSE)

Arguments

X.obs

Design matrix with missingness \(N \times p\)

y

Response vector \(N \times 1\)

pos_var

Index of selected covariates. The default is pos_var = 1:ncol(X.obs).

maxruns

Maximum number of iterations. The default is maxruns = 500.

tol_em

The tolerance to stop SAEM. The default is tol_em = 1e-7.

nmcmc

The MCMC length. The default is nmcmc = 2.

tau

Rate \(\tau\) in the step size \((k-k_{1})^{-\tau}\). The default is tau = 1.

k1

Number of first iterations \(k_{1}\) in the step size \((k-k_{1})^{-\tau}\). The default is k1=50.

seed

An integer as a seed set for the radom generator. The default value is 200.

print_iter

If TRUE, miss.saem will print the estimated parameters in each iteration of SAEM.

var_cal

If TRUE, miss.saem will calculate the variance of estimated parameters.

ll_obs_cal

If TRUE, miss.saem will calculate the observed log-likelihood.

Value

A list with components

mu

Estimated \(\mu\).

sig2

Estimated \(\Sigma\).

beta

Estiamated \(\beta\).

time_run

Execution time.

seqbeta

Sequence of \(\beta\) estimated in each iteration.

seqbeta_avg

Sequence of \(\beta\) with averaging in each iteration.

ll

Observed log-likelihood.

var_obs

Estimated variance for estimated parameters.

std_obs

Estimated standard error for estimated parameters.

Examples

# NOT RUN {
# Generate dataset
N <- 100  # number of subjects
p <- 3     # number of explanatory variables
mu.star <- rep(0,p)  # mean of the explanatory variables
Sigma.star <- diag(rep(1,p)) # covariance
beta.star <- c(1, 1,  0) # coefficients
beta0.star <- 0 # intercept
beta.true = c(beta0.star,beta.star)
X.complete <- matrix(rnorm(N*p), nrow=N)%*%chol(Sigma.star) +
              matrix(rep(mu.star,N), nrow=N, byrow = TRUE)
p1 <- 1/(1+exp(-X.complete%*%beta.star-beta0.star))
y <- as.numeric(runif(N)<p1)
# Generate missingness
p.miss <- 0.10
patterns <- runif(N*p)<p.miss #missing completely at random
X.obs <- X.complete
X.obs[patterns] <- NA

# SAEM
list.saem = miss.saem(X.obs,y)
print(list.saem$beta)
# }