tests/testthat/test-gr.R

In rlme: Rank-Based Estimation and Prediction in Random Effects Nested Models

context("GR method")

data('schools')

##### 3 Level GR #####

rlme.fit = rlme(y ~ 1 + sex + age + (1 | region) + (1 | region:school), schools, method="gr", rprpair='hl-disp')

test_that("GR calculates correct fixed effects for schools data", {  
  expect_equal(rlme.fit$fixed.effects$Estimate, c(0.1586624, -0.2953611, 0.2260327), tolerance=0.25)
})


test_that("GR calculates correct covariance matrix for school data", {
  correct.var.b = matrix(c(7.219013e-02, -0.0002000644, -1.672882e-05,
                           -2.000644e-04,  0.0114747092,  9.802962e-04,
                           -1.672882e-05,  0.0009802962,  2.629616e-02), nrow=3, ncol=3)
  
  expect_equal(unname(rlme.fit$var.b), correct.var.b, tolerance=0.25)
})


##### 3 Level HL-Disp #####

test_that("GR calculates correct random effect variances for school data (hl-disp)", {
  expect_equal(rlme.fit$random.effects$Variance, c(0.14703510, 0.01497416, 0.81229310), tolerance=0.25)
})


##### 3 Level Med-Mad #####

rlme.fit = rlme(y ~ 1 + sex + age + (1 | region) + (1 | region:school), schools, method="gr", rprpair='hl-disp')

test_that("GR calculates correct random effect variances for school data (med-mad)", {
  expect_equal(rlme.fit$random.effects$Variance, c(0.14703510, 0.01497416, 0.81229310), tolerance=0.25)
})


##### 2 Level GR #####

rlme.fit = rlme(y ~ 1 + sex + age + (1 | region), schools, method="gr", rprpair='hl-disp')

test_that("2 Level GR calculates correct fixed effects for schools data", {  
  # This is the expected value from the original package
  #expect_equal(rlme.fit$fixed.effects$Estimate, c(-0.09298845, -0.35939453, 0.11882249), tolerance=0.25)
  
  # when wilstep2 was switched over to minimize_dispersion, the values changed slightly
  # as the minimize_dispersion estimate was pretty close to the wwest estimate, I think the small
  # error is propogating through. 
  # I don't think it's a huge problem -- the test has been updated with the new value.
  
  expect_equal(rlme.fit$fixed.effects$Estimate, c(-0.08220573, -0.35643090, 0.12355131), tolerance=0.25)
})

test_that("2 Level GR calculates correct covariance matrix for school data", {
  # Original value (see comment in fixed effects test)
  #correct.var.b = matrix(c( 4.867082e-02, -0.0002499183, -5.661817e-05,
  #                         -2.499183e-04,  0.0080592483,  1.031766e-03,
  #                         -5.661817e-05,  0.0010317659,  1.867429e-02), nrow=3, ncol=3)
  
  correct.var.b = matrix(c( 4.853862e-02, -0.0002525695, -5.824018e-05,
                           -2.525695e-04,  0.0082960888,  1.067463e-03,
                           -5.824018e-05,  0.0010674635,  1.921460e-02), nrow=3, ncol=3)
  
  expect_equal(unname(rlme.fit$var.b), correct.var.b, tolerance = 0.01)
})


##### 2 Level HL-Disp #####

test_that("2 Level GR calculates correct random effect variances for school data (hl-disp)", {
  expect_equal(rlme.fit$random.effects$Variance, c(0.1409576, 0.8295030), tolerance=0.25)
})


##### 2 Level Med-Mad #####

rlme.fit = rlme(y ~ 1 + sex + age + (1 | region), schools, method="gr", rprpair='med-mad')

test_that("2 Level GR calculates correct random effect variances for school data (med-mad)", {
  # Value from original package
  # expect_equal(rlme.fit$random.effects$Variance, c(0.009342991, 0.718898005), tolerance=0.25)
  
  expect_equal(rlme.fit$random.effects$Variance, c(0.009548361, 0.726588754), tolerance=0.25)
})