rqr: Randomized Quantile Residuals for cub models.

In fhernanb/cubm: Tools to fit cub models via Maximum Likelihood

Description

This function generates Randomized Quantile Residuals.

Usage

rqr(y, pi, xi, m)

Arguments

y	vector with the response variable.
pi	value or vector with the values for pi parameter.
xi	value or vector with the values for pi parameter.
m	the maximum value.

Examples

# Example 1, without covariates -------------------------------------------

set.seed(456)
y <- rcub(n=1000, pi=0.15, xi=0.6, m=8)

mod1 <- cub(pi.fo=y~1, xi.fo=~1, m=8)
summary(mod1)

# Residuals for the correct model
r1 <- rqr(y, pi=mod1$fitted.pi, xi=mod1$fitted.xi, m=8)

# Residuals for a wrong model, using arbitrary estimated parameters
r2 <- rqr(y, pi=0.8, xi=0.2, m=8)

par(mfrow=c(1, 2))
car::qqPlot(r1, dist="norm", mean=0, sd=1, main='Correct model', las=1)
car::qqPlot(r2, dist="norm", mean=0, sd=1, main='Wrong model', las=1)

ks.test(r1, "pnorm", mean=0, sd=1)
ks.test(r2, "pnorm", mean=0, sd=1)


# Example 2, with covariates ----------------------------------------------

# Function to generate some random values from a CUB model
# pi explained by x1
# xi explained by x2

gendata <- function(n, m) {
  x1 <- runif(n)
  x2 <- rpois(n, lambda=4)
  pi <- pnorm(-1 + 2 * x1)
  xi <- pnorm( 4 - 1 * x2)
  y <- rcub(n=n, pi=pi, xi=xi, m=m)
  data.frame(y, x1, x2)
}

set.seed(12345)
dataset <- gendata(n=500, m=5)

# This is the correct model, pi explained by x1 and xi explained by x2
mod3 <- cub(pi.fo=y ~ x1, xi.fo= ~ x2, m=5, data=dataset)

# This is the wrong model
mod4 <- cub(pi.fo=y ~ x2, xi.fo=~ x1, m=5, data=dataset)

# The residuals for both models
r3 <- rqr(dataset$y, pi=mod3$fitted.pi, xi=mod3$fitted.xi, m=5)
r4 <- rqr(dataset$y, pi=mod4$fitted.pi, xi=mod4$fitted.xi, m=5)

par(mfrow=c(1, 2))
car::qqPlot(r3, dist="norm", mean=0, sd=1, main='Correct model', las=1)
car::qqPlot(r4, dist="norm", mean=0, sd=1, main='Wrong model', las=1)

ks.test(r3, "pnorm", mean=0, sd=1)
ks.test(r4, "pnorm", mean=0, sd=1)

fhernanb/cubm documentation built on Dec. 10, 2020, 1:24 p.m.