# tests/testthat/test-all.R In MNP: Fitting the Multinomial Probit Model

```rm(list=ls())
library(MNP)
library(testthat)
context("tests MNP")

test_that("tests MNP on the detergent data", {
# set random seed
set.seed(12345)

data(detergent)
# run the standard multinomial probit model with intercepts and the price
res1 <- mnp(choice ~ 1, choiceX = list(Surf=SurfPrice, Tide=TidePrice,
Wisk=WiskPrice, EraPlus=EraPlusPrice, Solo=SoloPrice, All=AllPrice),
cXnames = "price", data = detergent, n.draws = 30, burnin = 10,
thin = 3, verbose = TRUE)
# summarize the results
x <- summary(res1)
expect_that(length(x), is_equivalent_to(8))
expect_true("coef.table" %in% names(x))
expect_equal(x\$coef.table[4, 1], 1.864768, tolerance = 0.001)
expect_equal(x\$coef.table["(Intercept):Solo", "2.5%"], 1.108233, tolerance = 0.001)
expect_equal(x\$cov.table[10, 1], 0.9855489, tolerance = 0.001)
expect_equal(x\$cov.table["Tide:Tide", "mean"], 0.7952514, tolerance = 0.001)

# calculate the quantities of interest for the first 3 observations
x <- predict(res1, newdata = detergent[1:3,])
expect_that(length(x), is_equivalent_to(4))
expect_true("p" %in% names(x))
expect_that(dim(x\$o), is_equivalent_to(c(3, 6, 5)))
expect_equal(as.numeric(x\$p[1, 3]), 0.2, tolerance = 0.05)
expect_equal(as.numeric(x\$p[3, "Wisk"]), 0.4, tolerance = 0.05)
})

test_that("tests MNP on the Japanese election census", {
# set random seed
set.seed(12345)

# load the Japanese election data
data(japan)
# run the multinomial probit model with ordered preferences
res2 <- mnp(cbind(LDP, NFP, SKG, JCP) ~ gender + education + age, data = japan, verbose = TRUE)
# summarize the results
x <- summary(res2)
expect_that(length(x), is_equivalent_to(8))
expect_true("coef.table" %in% names(x))
expect_equal(x\$cov.table[2,1], 1.01725, tolerance = 0.01)

# calculate the predicted probabilities for the 10th observation
# averaging over 100 additional Monte Carlo draws given each of MCMC draw.
x <- predict(res2, newdata = japan[10,], type = "prob", n.draws = 100, verbose = TRUE)
expect_that(length(x), is_equivalent_to(2))
expect_true("p" %in% names(x))
expect_that(dim(x\$p), is_equivalent_to(c(1, 4, 5000)))
expect_that(x\$x[1, "age:LDP"], is_equivalent_to(50))
expect_that(x\$x[1, 1], is_equivalent_to(1))
})
```

## Try the MNP package in your browser

Any scripts or data that you put into this service are public.

MNP documentation built on June 22, 2024, 10:50 a.m.