Intro

Packages to install:

install.packages(c("lmtest", "ggplot2"))

Remarks on bptest()

Let’s make some data that is very heteroscedastic:

set.seed(4040)
observations <- 500
trumpet <-
  data.frame(
    x = 1:observations,
    y = rnorm(observations, mean = 0, sd = 1:observations)
  )

You can see the heteroscedasticity:

ggplot(trumpet, aes(x = x, y = y)) +
  geom_point()

Now have lmtest::bptest() do its work:

mod <- lm(y ~ x, data = trumpet)
bptest(mod)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod
## BP = 100.15, df = 1, p-value < 2.2e-16

But the above P-value is only human-readable representation of the P-value that the routine computed. To see that this is so, look into the structure of the object produced by a call to bptest():

results <- bptest(mod)
str(results)

## List of 5
##  $ statistic: Named num 100
##   ..- attr(*, "names")= chr "BP"
##  $ parameter: Named num 1
##   ..- attr(*, "names")= chr "df"
##  $ method   : chr "studentized Breusch-Pagan test"
##  $ p.value  : Named num 1.41e-23
##   ..- attr(*, "names")= chr "BP"
##  $ data.name: chr "mod"
##  - attr(*, "class")= chr "htest"

Aha! We can access the P-value as follows:

results$p.value

##           BP 
## 1.413889e-23

So I figure that the P-value you computed by hand is probably the same as what bptest() got, under the hood.

Choosing R Packages

Any open-source software of any note accumulates an overgrown jungle of packages. People do attempt to curate lists of the best packages for given tasks. One highly-regarded (but also very incomplete) set of sets of packages is provided in CRAN’s Task Views.