Running examples in 'qrcm-Ex.R' failed
The error most likely occurred in:

> ### Name: iqr
> ### Title: Quantile Regression Coefficients Modeling
> ### Aliases: iqr
> ### Keywords: models regression
> 
> ### ** Examples
> 
> 
>   ##### Using simulated data in all examples
> 
> 
>   ##### Example 1
>   
>   n <- 1000
>   x <- runif(n)
>   y <- rnorm(n, 1 + x, 1 + x)
>   # true quantile function: Q(p | x) = beta0(p) + beta1(p)*x, with 
>     # beta0(p) = beta1(p) = 1 + qnorm(p)
>                               
>   # fit the true model: b(p) = (1 , qnorm(p))
>   m1 <- iqr(y ~ x, formula.p = ~ I(qnorm(p)))
>   # the fitted quantile regression coefficient functions are
>     # beta0(p) = m1$coef[1,1] + m1$coef[1,2]*qnorm(p)
>     # beta1(p) = m1$coef[2,1] + m1$coef[2,2]*qnorm(p)
> 
>   # a basis b(p) = (1, p), i.e., beta(p) is assumed to be a linear function of p
>   m2 <- iqr(y ~ x, formula.p = ~ p)
> 
>   # a 'rich' basis b(p) = (1, p, p^2, log(p), log(1 - p))
>   m3 <- iqr(y ~ x, formula.p = ~ p + I(p^2) + I(log(p)) + I(log(1 - p)))
> 
>   # 'slp' creates an orthogonal spline basis using shifted Legendre polynomials
>   m4 <- iqr(y ~ x, formula.p = ~ slp(p, k = 3)) # note that this is the default
>   
>   # 'plf' creates the basis of a piecewise linear function
>   m5 <- iqr(y ~ x, formula.p = ~ plf(p, knots = c(0.1,0.9)))
> 
>   
>   summary(m1)

Call: iqr(formula = y ~ x, formula.p = ~I(qnorm(p)))

converged: TRUE 
n. of iterations: 11 
n. of observations: 1000 
n. of free parameters: 4 

###################### 
###################### 

Coefficients:
             (Intercept)  I(qnorm(p))
(Intercept)  0.9737       1.0630     
x            0.9707       0.9837     

Standard errors:
             (Intercept)  I(qnorm(p))
(Intercept)  0.08412      0.06137    
x            0.17485      0.13048    

###################### 
###################### 

Wald test for x:
            chi-square df P(> chi)    
(Intercept)     408.09  2   <2e-16 ***
x                81.45  2   <2e-16 ***


Wald test for b(p):
            chi-square df P(> chi)    
(Intercept)      836.6  2   <2e-16 ***
I(qnorm(p))     1680.9  2   <2e-16 ***

Minimized loss function: 436.9049>   summary(m1, p = c(0.25,0.5,0.75))

Call: iqr(formula = y ~ x, formula.p = ~I(qnorm(p)))

p = 0.25 

Coefficients:
            Estimate std.err z value p(>|z|))   
(Intercept)  0.25675 0.09113   2.817  0.00484 **
x            0.30724 0.18931   1.623  0.10461   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



p = 0.5 

Coefficients:
            Estimate std.err z value p(>|z|))    
(Intercept)  0.97372 0.08412  11.576  < 2e-16 ***
x            0.97073 0.17485   5.552 2.83e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



p = 0.75 

Coefficients:
            Estimate std.err z value p(>|z|))    
(Intercept)   1.6907  0.0963  17.557  < 2e-16 ***
x             1.6342  0.2020   8.091 5.91e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



>   test.fit(m1)
                    statistic p-value
Kolmogorov-Smirnov 0.01613508    0.68
Cramer-Von Mises   0.04139944    0.55
>   par(mfrow = c(1,2)); plot(m1, ask = FALSE)
>   # see the documentation for 'summary.iqr', 'test.fit', and 'plot.iqr'
> 
> 
> 
> 
>   ##### Example 2 ### excluding coefficients
>   
>   n <- 1000
>   x <- runif(n)
>   qy <- function(p,x){(1 + qnorm(p)) + (1 + log(p))*x}
>   # true quantile function: Q(p | x) = beta0(p) + beta1(p)*x, with
>     # beta0(p) = 1 + qnorm(p) 
>     # beta1(p) = 1 + log(p)
>   
>   y <- qy(runif(n), x) # to generate y, plug uniform p in qy(p,x) 
>   iqr(y ~ x, formula.p = ~ I(qnorm(p)) + I(log(p)))

Call:
iqr(formula = y ~ x, formula.p = ~I(qnorm(p)) + I(log(p)))

Coefficients:
             (Intercept)  I(qnorm(p))  I(log(p))
(Intercept)  0.98337      0.97909      0.01704  
x            0.81548      0.27308      0.79011  

> 
>   # I would like to exclude log(p) from beta0(p), and qnorm(p) from beta1(p)
>   # I set to 0 the corresponding entries of 's'
> 
>   s <- matrix(1,2,3); s[1,3] <- s[2,2] <- 0
>   iqr(y ~ x, formula.p = ~ I(qnorm(p)) + I(log(p)), s = s)

Call:
iqr(formula = y ~ x, formula.p = ~I(qnorm(p)) + I(log(p)), s = s)

Coefficients:
             (Intercept)  I(qnorm(p))  I(log(p))
(Intercept)  0.9683       1.0297       0.0000   
x            1.0244       0.0000       1.0353   

> 
> 
> 
> 
>   ##### Example 3 ### excluding coefficients when b(p) is singular
>   
>   n <- 1000
>   x <- runif(n)
>   qy <- function(p,x){(1 + log(p) - 2*log(1 - p)) + (1 + log(p/(1 - p)))*x} 
>   # true quantile function: Q(p | x) = beta0(p) + beta1(p)*x, with
>     # beta0(p) = 1 + log(p) - 2*log(1 - p)
>     # beta1(p) = 1 + log(p/(1 - p))
> 
>   y <- qy(runif(n), x) # to generate y, plug uniform p in qy(p,x)
> 
>   iqr(y ~ x, formula.p = ~ I(log(p)) + I(log(1 - p)) + I(log(p/(1 - p))))

Call:
iqr(formula = y ~ x, formula.p = ~I(log(p)) + I(log(1 - p)) + 
    I(log(p/(1 - p))))

Coefficients:
             (Intercept)  I(log(p))  I(log(1 - p))  I(log(p/(1 - p)))
(Intercept)   0.8254       0.9094    -1.7568         0.0000          
x             1.2195       1.0247    -1.2111         0.0000          

>   # log(p/(1 - p)) is dropped due to singularity
>   
>   # I want beta0(p) to be a function of log(p) and log(1 - p),
>   # and beta1(p) to depend on log(p/(1 - p)) alone
> 
>   s <- matrix(1,2,4); s[2,2:3] <- 0
>   iqr(y ~ x, formula.p = ~ I(log(p)) + I(log(1 - p)) + I(log(p/(1 - p))), s = s)

Call:
iqr(formula = y ~ x, formula.p = ~I(log(p)) + I(log(1 - p)) + 
    I(log(p/(1 - p))), s = s)

Coefficients:
             (Intercept)  I(log(p))  I(log(1 - p))  I(log(p/(1 - p)))
(Intercept)   0.7548       0.8689    -1.8002         0.0000          
x             1.3770       0.0000     0.0000         1.1151          

>   # log(p/(1 - p)) is not dropped
> 
> 
> 
> 
>   ##### Example 4 ### using slp to test deviations from normality
>   
>   n <- 1000
>   x <- runif(n)
>   y <- rnorm(n, 2 + x) 
>   # the true model is normal, i.e., b(p) = (1, qnorm(p))
>   
>   summary(iqr(y ~ x, formula.p = ~ I(qnorm(p)) + slp(p,3))) 

Call: iqr(formula = y ~ x, formula.p = ~I(qnorm(p)) + slp(p, 3))

converged: TRUE 
n. of iterations: 8 
n. of observations: 1000 
n. of free parameters: 10 

###################### 
###################### 

Coefficients:
             (Intercept)  I(qnorm(p))  slp(p, 3)slp1  slp(p, 3)slp2
(Intercept)   2.07036      1.01906     -0.12741        0.10301     
x             0.51578     -0.28156      0.49475       -0.06645     
             slp(p, 3)slp3
(Intercept)   0.08096     
x             0.01457     

Standard errors:
             (Intercept)  I(qnorm(p))  slp(p, 3)slp1  slp(p, 3)slp2
(Intercept)  0.6040       0.2892       0.4877         0.0695       
x            1.0482       0.5028       0.8471         0.1208       
             slp(p, 3)slp3
(Intercept)  0.1257       
x            0.2197       

###################### 
###################### 

Wald test for x:
            chi-square df P(> chi)    
(Intercept)     1422.8  5   <2e-16 ***
x                102.4  5   <2e-16 ***


Wald test for b(p):
            chi-square df P(> chi)    
(Intercept)     53.321  2 2.64e-12 ***
I(qnorm(p))     33.223  2 6.11e-08 ***
slp(p, 3)        9.701  6    0.138    

Minimized loss function: 268.5465>   # if slp(p,3) is not significant, no deviation from normality
> 
> 
> 
> 
>   ##### Example 5 ### formula without intercept
>   
>   n <- 1000
>   x <- runif(n)
>   y <- runif(n, 0,x) 
> 
>   # True quantile function: Q(p | x) = p*x, i.e., beta0(p) = 0, beta1(p) = p
>   # When x = 0, all quantiles of y are 0, i.e., the distribution is degenerated
>   # To explicitly model this, remove the intercept from 'formula'
>   
>   iqr(y ~ -1 + x, formula.p = ~ p)

Call:
iqr(formula = y ~ -1 + x, formula.p = ~p)

Coefficients:
   (Intercept)  p       
x  -0.00197      1.00820

>   
>   # the true model does not have intercept in b(p) either:
> 
>   iqr(y ~ -1 + x, formula.p = ~ -1 + p)

Call:
iqr(formula = y ~ -1 + x, formula.p = ~-1 + p)

Coefficients:
   p    
x  1.005

> 
> 
> 
> 
>   ##### Example 6 ### no covariates, strictly positive outcome
>   
>   n <- 1000
>   y <- rgamma(n, 3,1) 
> 
>   # you know that Q(0) = 0
>   # remove intercept from 'formula.p', and use b(p) such that b(0) = 0
>   
>   summary(iqr(y ~ 1, formula.p = ~ -1 + slp(p,5))) # shifted Legendre polynomials

Call: iqr(formula = y ~ 1, formula.p = ~-1 + slp(p, 5))

converged: TRUE 
n. of iterations: 6 
n. of observations: 1000 
n. of free parameters: 5 

###################### 
###################### 

Coefficients:
             slp1    slp2    slp3    slp4    slp5  
(Intercept)  2.6823  0.6965  0.8190  0.2034  0.3222

Standard errors:
             slp1     slp2     slp3     slp4     slp5   
(Intercept)  0.07201  0.06873  0.06248  0.05851  0.04163

###################### 
###################### 

Wald test for b(p):
     chi-square df P(> chi)    
slp1    1387.54  1  < 2e-16 ***
slp2     102.71  1  < 2e-16 ***
slp3     171.82  1  < 2e-16 ***
slp4      12.08  1 0.000509 ***
slp5      59.90  1 9.99e-15 ***

Minimized loss function: 447.0672>   summary(iqr(y ~ 1, formula.p = ~ -1 + sin(p*pi/2) + I(qbeta(p,2,4)))) # unusual basis

Call: iqr(formula = y ~ 1, formula.p = ~-1 + sin(p * pi/2) + I(qbeta(p, 
    2, 4)))

converged: TRUE 
n. of iterations: 9 
n. of observations: 1000 
n. of free parameters: 2 

###################### 
###################### 

Coefficients:
             sin(p * pi/2)  I(qbeta(p, 2, 4))
(Intercept)  -1.332         11.327           

Standard errors:
             sin(p * pi/2)  I(qbeta(p, 2, 4))
(Intercept)  0.2922         0.6039           

###################### 
###################### 

Wald test for b(p):
                  chi-square df P(> chi)    
sin(p * pi/2)          20.76  1  5.2e-06 ***
I(qbeta(p, 2, 4))     351.86  1  < 2e-16 ***

Minimized loss function: 447.1758>   summary(iqr(y ~ 1, formula.p = ~ -1 + I(sqrt(p))*I(log(1 - p)))) # you can include interactions

Call: iqr(formula = y ~ 1, formula.p = ~-1 + I(sqrt(p)) * I(log(1 - 
    p)))

converged: TRUE 
n. of iterations: 9 
n. of observations: 1000 
n. of free parameters: 3 

###################### 
###################### 

Coefficients:
             I(sqrt(p))  I(log(1 - p))  I(sqrt(p)):I(log(1 - p))
(Intercept)   4.378       3.927         -4.612                  

Standard errors:
             I(sqrt(p))  I(log(1 - p))  I(sqrt(p)):I(log(1 - p))
(Intercept)  0.3867      1.1660         1.0815                  

###################### 
###################### 

Wald test for b(p):
                         chi-square df P(> chi)    
I(sqrt(p))                   128.17  1  < 2e-16 ***
I(log(1 - p))                 11.34  1 0.000759 ***
I(sqrt(p)):I(log(1 - p))      18.18  1 2.01e-05 ***

Minimized loss function: 446.9422> 
> 
> 
> 
>   ##### Example 7 ### revisiting the classical linear model
>   
>   n <- 1000
>   x <- runif(n)
>   y <- 2 + 3*x + rnorm(n,0,1) # beta0 = 2, beta1 = 3
>   
>   iqr(y ~ x, formula.p = ~ I(qnorm(p)), s = matrix(c(1,1,1,0),2))

Call:
iqr(formula = y ~ x, formula.p = ~I(qnorm(p)), s = matrix(c(1, 
    1, 1, 0), 2))

Coefficients:
             (Intercept)  I(qnorm(p))
(Intercept)  2.1253       0.9959     
x            2.8947       0.0000     

>   # first column of coefficients: (beta0, beta1)
>   # top-right coefficient: residual standard deviation
>   
> 
> 
> 
>   ##### Example 8 ### censored data
>   
>   n <- 1000
>   x <- runif(n,0,5)
> 	
>   u <- runif(n)
>   beta0 <- -log(1 - u)
>   beta1 <- 0.2*log(1 - u)
>   t <- beta0 + beta1*x  # time variable
>   c <- rexp(n,2)        # censoring variable
>   y <- pmin(t,c)        # observed events
>   d <- (t <= c)         # 1 = event, 0 = censored
>   
>   iqr(Surv(y,d) ~ x, formula.p = ~ I(log(1 - p)))

Call:
iqr(formula = Surv(y, d) ~ x, formula.p = ~I(log(1 - p)))

Coefficients:
             (Intercept)  I(log(1 - p))
(Intercept)   0.008643    -0.871372    
x            -0.001681     0.172578    

>   
>   ##### Example 8 (cont.) ### censored and truncated data
> 
>   z <- rexp(n,10)   # truncation variable
>   w <- which(y > z) # only observe z,y,d,x when y > z
>   z <- z[w]; y <- y[w]; d <- d[w]; x <- x[w]
> 
>   iqr(Surv(z,y,d) ~ x, formula.p = ~ I(log(1 - p)))
Error in aeqSurv(Y) : 
  aeqSurv exception, an interval has effective length 0
Calls: iqr ... ctiqr.internal -> km -> survfit -> coxph -> aeqSurv
Execution halted