context("Testing blm fit correctness")
seed <- as.integer(1000 * rnorm(1))
test_that(paste("Deviance of a blm fit can max be equal to lm fit using seed", seed), {
w0 <- 0.3 ; w1 <- 1.1 ; b <- 1.3
x <- rnorm(50)
y <- rnorm(50, w1 * x + w0, 1/b)
mod <- blm(y~x, beta=b, data=data.frame(x=x, y=y))
expect_gte(deviance(mod), deviance(lm(y~x, data=data.frame(x=x, y=y))))
})
seed <- as.integer(1000 * rnorm(1))
test_that(paste("1/var of residuals close to precision of data using seed", seed), {
w0 <- 0.3 ; w1 <- 1.1 ; b <- 0.1
x <- rnorm(1000)
y <- rnorm(1000, w1 * x + w0, 1/sqrt(b))
mod1000 <- blm(y~x, data=data.frame(x=x, y=y))
expect_equal(1/var(resid(mod1000)), b, tolerance=.1)
})
seed <- as.integer(1000 * rnorm(1))
test_that(paste("Covariance of coefficients of a blm fit reflects the precision of the model using seed", seed), {
w0 <- 0.3 ; w1 <- 1.1 ; b <- 1.3
x <- rnorm(100)
y <- rnorm(100, w1 * x + w0, 1/b)
mod13 <- blm(y~x, beta=b, data=data.frame(x=x, y=y))
covar.blm(mod13)
b33 <- 3.3
y33 <- rnorm(100, w1 * x + w0, 1/b33)
mod33 <- blm(y33~x, beta=b33, data=data.frame(x=x, y=y33))
covar.blm(mod33)
#model with more precision should have smaller variance on the coefficients
expect_equal(all(as.vector(abs(covar.blm(mod13))) >
as.vector(abs(covar.blm(mod33)))),
TRUE)
#model provided with higher precision should have less variance on the
#coefficients than equal model provided with higher precision
mod13b <- blm(y~x, beta=b33, data=data.frame(x=x, y=y))
expect_equal(all(as.vector(abs(covar.blm(mod13))) >
as.vector(abs(covar.blm(mod13b)))),
TRUE)
})
seed <- as.integer(1000 * rnorm(1))
test_that(paste("Variance of predicted values increases with distance to data points using seed", seed), {
w0 <- 0.3 ; w1 <- 1.1 ; b <- 1.3
x <- rnorm(100)
y <- rnorm(100, w1 * x + w0, 1/b)
mod <- blm(y~x, beta=b, data=data.frame(x=x, y=y))
x2 <- 0:100
y2 <- predict(mod, data.frame(x=x2), se.fit=T)
expect_equal(cor(x2, y2$se.fit, method='spearman'),1)
})
seed <- as.integer(1000 * rnorm(1))
test_that(paste("Roughly 66% of random sampled data using model parameters falls into fitted values +/- sqrt(var) using seed", seed), {
w0 <- 1
w1 <- 2
x <- seq(-100,100,10)
b <- 0.001
y <- w0 + w1*x + rnorm(length(x), mean=0, sd=sqrt(1/b) )
blm_mod_with_b <- blm(y~x, beta = b)
test <- function(){
y <- w0 + w1*x + rnorm(length(x), mean=0, sd=sqrt(1/b) )
fits <- fitted(blm_mod_with_b, var=T)
fits$sd <- sqrt(fits$var)
pos <- vapply(seq_along(y), function(i) ifelse( (y[i] >= fits$mean[i] - fits$sd[i]) &
(y[i] <= fits$mean[i] + fits$sd[i]),
TRUE, FALSE), TRUE )
sum(pos)/length(pos)
}
tests <- replicate(10000, test())
expect_equal(mean(tests), .65, tolerance = .1)
})
seed <- as.integer(1000 * rnorm(1))
test_that(paste("Self-updating using same data improves fit using seed", seed), {
w0 <- 0.3 ; w1 <- 1.1 ; b <- 1.3
x <- rnorm(100)
y <- rnorm(100, w1 * x + w0, 1/b)
mod1 <- blm(y~x)
mod2 <- update(mod1)
mod3 <- update(mod2)
mod4 <- update(mod3)
dist_4_vs_4 <- mahal(mod4$posterior, coef(mod4))
dist_3_vs_4 <- mahal(mod4$posterior, coef(mod3))
dist_2_vs_4 <- mahal(mod4$posterior, coef(mod2))
dist_1_vs_4 <- mahal(mod4$posterior, coef(mod1))
dist_prior_vs_4 <- mahal(mod4$posterior, mod1$prior$means)
expect_equal(dist_4_vs_4, 0)
expect_lte(dist_3_vs_4, dist_2_vs_4)
expect_lte(dist_2_vs_4, dist_1_vs_4)
expect_lte(dist_1_vs_4, dist_prior_vs_4)
dev4 <- deviance(mod4)
dev3 <- deviance(mod3)
dev2 <- deviance(mod2)
dev1 <- deviance(mod1)
expect_lte(dev4, dev3)
expect_lte(dev3, dev2)
expect_lte(dev2, dev1)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.