Nothing
context("random_effects")
# creating cluster variable for test purpose only
set.seed(1)
# Weibull parameters
shape <- 0.9
scale <- 2
# number of simulated datasets
NFile <- 1
# number of individuals per dataset
n <- 2000
# number of clusters
NCluster <- 20
# data frame
data.rd <- data.frame(cluster=seq(1:NCluster))
cluster <- sample(rep(1:NCluster,each=n/NCluster))
don <- data.frame(num=1:n, cluster=factor(cluster)) # be careful, cluster needs to be a factor !
don <- merge(don, data.rd, by="cluster")[, union(names(don), names(data.rd))]
don <- don[order(don$num),]
rownames(don) <- NULL
# theoretical standard deviation
sd1 <- 0.1
# maximum follow-up time
max.time <- 5
wj <- rnorm(NCluster,mean=0,sd=sd1)
don$wj <- wj[don$cluster]
# simulated times
u <- runif(n)
don$fu <- exp( 1/shape*(log(-log(1-u)) - don$wj) + log(scale))
# censoring
don$dead <- ifelse(don$fu <= max.time,1,0)
don$fu <- pmin(don$fu,max.time)
# fitting
mod.rd <- survPen(~smf(fu)+rd(cluster),data=don,t1=fu,event=dead,detail.beta=TRUE,detail.rho=TRUE)
test_that("random effects models work", {
expect_true(inherits(mod.rd,"survPen"))
})
test_that("estimated variance is ok", {
expect_true(abs(exp(summary(mod.rd)$random.effects)[1]
- exp(-0.5*log(mod.rd$lambda)-0.5*log(mod.rd$S.scale))[2]) < 1e-08)
})
test_that("prediction ok", {
expect_false(predict(mod.rd,data.frame(fu=5,cluster=1))$surv == predict(mod.rd,data.frame(fu=5,cluster=1),exclude.random=TRUE)$surv)
})
test_that("penalized likelihoods from model and summary matches", {
expect_equal(summary(mod.rd)$penalized.likelihood , mod.rd$ll.pen)
})
test_that("summary has right dimensions", {
expect_true(all(dim(summary(mod.rd)$random.effects)==c(1,2)))
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.