R Under development (unstable) (2024-02-01 r85851 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> require("expectreg")
Loading required package: expectreg
Loading required package: parallel
Loading required package: mboost
Loading required package: stabs
Loading required package: BayesX
Loading required package: shapefiles
Loading required package: foreign

Attaching package: 'shapefiles'

The following objects are masked from 'package:foreign':

    read.dbf, write.dbf

Loading required package: Matrix
> 
> set.seed(9484470)
> 
> ###compare 0.5-expectile regression with lm in appropiate cases
> 
> ma <- matrix(runif(200)*10,nrow=100)
> 
> testfunction <- function(mod1,mod2) {
+     stopifnot(max(abs(fitted(mod1)-fitted(mod2)))<0.5)
+ }
> 
> #univariate linear model
> dat <- data.frame(regressand = 3*ma[,1],covariate.1=ma[,2])
> m1 <- expectreg.ls(regressand~ covariate.1,data=dat,smooth="f",expectiles=c(0.5))
> #m2 <- expectile.sheets(regressand~covariate.1,data=dat,smooth="f",expectiles=c(0.5))
> 
> m3 <- expectreg.ls(regressand~ covariate.1,data=dat,smooth="f",expectiles=c(0.5),estimate="restricted")
> m4 <- expectreg.ls(regressand~ covariate.1,data=dat,smooth="f",expectiles=c(0.5),estimate="bundle")
> 
> m5 <- expectreg.boost(regressand~ bols(covariate.1),data=dat,expectiles=c(0.5),mstop=112,cv=F)
Expectile  0.5 
> 
> mcomp <- lm(regressand~covariate.1,data=dat)
> 
> testfunction(m1,mcomp)
> #testfunction(m2,mcomp)
> testfunction(m3,mcomp)
> testfunction(m4,mcomp)
> testfunction(m5,mcomp)
> 
> # univariate model involving sin transformation
> 
> dat <- data.frame(regressand = sin(ma[,1]),covariate=ma)
> m1 <- expectreg.ls(regressand~rb( covariate.1,"pspline"),data=dat,smooth="f",expectiles=c(0.5),lambda=0.000000001)
> 
> #m2 <- expectile.sheets(regressand~rb( covariate.1,"pspline"),data=dat,smooth="acv",expectiles=c(0.5))
> 
> m3 <- expectreg.ls(regressand~rb( covariate.1,"pspline"),data=dat,smooth="schall",expectiles=c(0.5),estimate="restricted")
> m4 <- expectreg.ls(regressand~rb( covariate.1,"pspline"),data=dat,smooth="schall",expectiles=c(0.5),estimate="bundle")
> 
> m5 <- expectreg.boost(regressand~bbs( covariate.1),data=dat,expectiles=c(0.5),mstop=500,cv=F)
Expectile  0.5 
> 
> 
> mcomp <- lm(regressand~sin(covariate.1),data=dat)
> 
> testfunction(m1,mcomp)
> #testfunction(m2,mcomp)
> testfunction(m3,mcomp)
> testfunction(m4,mcomp)
> testfunction(m5,mcomp)
> 
> 
> # bivariate model: linear and sin
> dat <- data.frame(regressand = sin(ma[,1])+3*ma[,2],covariate=ma)
> 
> m1 <- expectreg.ls(regressand~rb( covariate.1,"pspline")+rb( covariate.2,"pspline"),data=dat,smooth="f",expectiles=c(0.5),lambda=0.000000001)
> #m2 <- expectile.sheets(regressand~rb( covariate.1,"pspline")+base( covariate.2,"pspline"),data=dat,smooth="acv",expectiles=c(0.5))
> 
> m3 <- expectreg.ls(regressand~rb( covariate.1,"pspline")+rb( covariate.2,"pspline"),data=dat,smooth="schall",expectiles=c(0.5),estimate="restricted")
> m4 <- expectreg.ls(regressand~rb( covariate.1,"pspline")+rb( covariate.2,"pspline"),data=dat,smooth="schall",expectiles=c(0.5),estimate="bundle")
> 
> m5 <- expectreg.boost(regressand~bbs( covariate.1)+bbs(covariate.2),data=dat,expectiles=c(0.5),mstop=1200,cv=F)
Expectile  0.5 
> 
> mcomp <- lm(regressand~sin(covariate.1)+covariate.2,data=dat)
> 
> testfunction(m1,mcomp)
> #testfunction(m2,mcomp)
> testfunction(m3,mcomp)
> testfunction(m4,mcomp)
> testfunction(m5,mcomp)
> 
> 
> set.seed(9484470)
> data(dutchboys)
> dutchb <- dutchboys[sample(nrow(dutchboys),150),]
> 
> expect <- c(0.05,0.5,0.95)
> 
> 
> ###Values
> 
> 
> mLaws <- expectreg.ls(hgt~rb(age,"pspline")+rb(wgt,"pspline"),data=dutchb,smooth="gcv",expectiles=expect)
> mSheets <- expectreg.ls(hgt~rb(age,"pspline")+rb(wgt,"pspline"),data=dutchb,smooth="f",lambda=1,expectiles=expect,estimate="sheets")
> mNoncross <-  expectreg.qp(hgt~rb(age,"pspline"),data=dutchb,smooth="f",expectiles=expect)
> mBundle <- expectreg.ls(hgt~rb(age,"pspline")+rb(wgt,"pspline"),data=dutchb,smooth="schall",expectiles=expect,estimate="bundle")
> mRestricted <- expectreg.ls(hgt~rb(age,"pspline")+rb(wgt,"pspline"),data=dutchb,smooth="schall",expectiles=expect,estimate="restricted")
> mBoost <- expectreg.boost(hgt~bbs(age)+bbs(wgt),data=dutchb,expectiles=expect,mstop=400)
Expectile  0.05 , Boosting iterations:  400 
Expectile  0.5 , Boosting iterations:  96 
Expectile  0.95 , Boosting iterations:  378 
> 
> ###INTERCEPTS
> stopifnot(length(mLaws$intercepts) == length(mLaws$asymmetries))
> #stopifnot(length(mSheets$intercepts) == length(mSheets$asymmetries))
> stopifnot(length(mNoncross$intercepts) == length(mNoncross$asymmetries))
> stopifnot(length(mBundle$intercepts) == length(mBundle$asymmetries))
> stopifnot(length(mRestricted$intercepts) == length(mRestricted$asymmetries))
> 
> ###COEFFICIENTS A matrix of all the coefficients, for each base element a row and for each expectile a column. 
> stopifnot(length(mLaws$coefficients)==2)
> #stopifnot(length(mSheets$coefficients)==2)
> stopifnot(length(mNoncross$coefficients)==1)
> stopifnot(length(mBundle$coefficients)==2)
> stopifnot(length(mRestricted$coefficients)==2)
> 
> for(i in 1:2) {
+     stopifnot(nrow(mLaws$coefficients[[i]])== 21)
+     stopifnot(ncol(mLaws$coefficients[[i]])== length(mLaws$asymmetries))
+ }
> #for(i in 1:2) {
> #   stopifnot(ncol(mSheets$coefficients[[i]])==length(mSheets$asymmetries))
> #}
> 
> 
> stopifnot(ncol(mNoncross$coefficients[[1]])== length(mNoncross$asymmetries))
> 
> 
> for(i in 1:2) {
+     stopifnot(nrow(mBundle$coefficients[[i]])== 21)
+     stopifnot(ncol(mBundle$coefficients[[i]])== length(mBundle$asymmetries))
+ }
> for(i in 1:2) {
+     stopifnot(nrow(mRestricted$coefficients[i])== 21)
+     stopifnot(ncol(mRestricted$coefficients[i])== length(mRestricted$asymmetries))
+ }
> 
> ###VALUES
> 
> ###RESPONSE
> stopifnot(isTRUE(all.equal(dutchb$hgt, mLaws$response,check.attributes=F)))
> stopifnot(isTRUE(all.equal(dutchb$hgt, mSheets$response,check.attributes=F)))
> stopifnot(isTRUE(all.equal(dutchb$hgt, mNoncross$response,check.attributes=F)))
> stopifnot(isTRUE(all.equal(dutchb$hgt, mBundle$response,check.attributes=F)))
> stopifnot(isTRUE(all.equal(dutchb$hgt, mRestricted$response,check.attributes=F)))
> stopifnot(isTRUE(all.equal(dutchb$hgt, mBoost$response,check.attributes=F)))
> 
> ###COVARIATES
> stopifnot(isTRUE(all.equal(dutchb$age, mLaws$covariates$age )))
> stopifnot(isTRUE(all.equal(dutchb$age, mSheets$covariates$age )))
> stopifnot(isTRUE(all.equal(dutchb$age, mNoncross$covariates$age )))
> stopifnot(isTRUE(all.equal(dutchb$age, mBundle$covariates$age)))
> stopifnot(isTRUE(all.equal(dutchb$age, mRestricted$covariates$age)))
> 
> ###EXPECTILES
> stopifnot(isTRUE(all.equal(mLaws$asymmetries, expect)))
> stopifnot(isTRUE(all.equal(mSheets$asymmetries, expect)))
> stopifnot(isTRUE(all.equal(mNoncross$asymmetries, expect)))
> stopifnot(isTRUE(all.equal(mBundle$asymmetries, expect)))
> stopifnot(isTRUE(all.equal(mRestricted$asymmetries, expect)))
> stopifnot(isTRUE(all.equal(mBoost$asymmetries, expect)))
> 
> ###EFFECTS
> stopifnot(length(mLaws$covariates)==length(mLaws$effects))
> stopifnot(length(mSheets$covariates)==length(mSheets$effects))
> stopifnot(length(mNoncross$covariates)==length(mNoncross$effects))
> stopifnot(length(mBundle$covariates)==length(mBundle$effects))
> stopifnot(length(mRestricted$covariates)==length(mRestricted$effects))
> stopifnot(length(mBoost$covariates)==length(mBoost$effects))
> 
> 
> #Methods
> 
> ###predict
> 
> stopifnot(isTRUE(all.equal((predict(mLaws,newdata = dutchb))$fitted, fitted(mLaws))))
> #stopifnot(isTRUE(all.equal((predict(mSheets,newdata = dutchb))$fitted, fitted(mSheets))))
> #stopifnot(isTRUE(all.equal((predict(mNoncross,newdata = dutchb))$fitted, fitted(mNoncross))))
> 
> stopifnot(isTRUE(all.equal((predict(mBundle,newdata = dutchb))$fitted, fitted(mBundle))))
> stopifnot(isTRUE(all.equal((predict(mRestricted,newdata = dutchb))$fitted, fitted(mRestricted))))
> 
> stopifnot(isTRUE(all.equal((predict(mBoost,newdata = dutchb))$fitted, fitted(mBoost))))
> 
> ###resid
> stopifnot(isTRUE(all.equal(resid(mLaws)$fitted, mLaws$response-mLaws$fitted)))
> stopifnot(isTRUE(all.equal(resid(mSheets)$fitted, mSheets$response-mSheets$fitted)))
> stopifnot(isTRUE(all.equal(resid(mNoncross)$fitted, mNoncross$response-mNoncross$fitted)))
> 
> stopifnot(isTRUE(all.equal(resid(mBundle)$fitted, mBundle$response-mBundle$fitted)))
> stopifnot(isTRUE(all.equal(resid(mRestricted)$fitted, mRestricted$response-mRestricted$fitted)))
> 
> stopifnot(isTRUE(all.equal(resid(mBoost)$fitted, mBoost$response-mBoost$fitted)))
> 
> proc.time()
   user  system elapsed 
   6.39    0.62    7.01