R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Linear Regression Based on Identity, Spatial Sign or Spatial...

mv.l1lm

R Documentation

Linear Regression Based on Identity, Spatial Sign or Spatial Rank Scores

Description

This function fits a multivariate linear regression model based on identity, spatial sign or spatial rank scores. Both inner and outer standardization are possible.

Usage

mv.l1lm(formula, scores = "identity", stand = "outer", 
        maxiter = 1000, eps = 1e-06, eps.S = 1e-06, 
        x = TRUE, y = TRUE, data, subset, na.action)

Arguments

`formula`	an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. The left part of the formula (the response) must be a n x p matrix with at least two columns.
`scores`	score to be used. Can be either `"identity"`, `"sign"` or `"rank"`.
`stand`	can be `"outer"` or `"inner"`.
`maxiter`	maximum number of iterations. Used only for `score = "sign"` and `score = "rank"`.
`eps`	convergence tolerance. Used only for `score = "sign"` or `score = "rank"`.
`eps.S`	lower limit for the residual norms. Used only for `score = "sign"` or `score = "rank"` in the iteration procedure to avoid to divide by a zero norm.
`x`	logical. Indicating whether the design matrix 'x' returned from the model matrix should be stored. Default is `TRUE`. Might be needed for example in the `anova` function.
`y`	logical. Indicating whether the response matrix 'y' should be stored. Default is `TRUE`.
`data`	an optional data frame, list or environment (or object coercible by 'as.data.frame' to a data frame) containing the variables in the model. If not found in 'data', the variables are taken from 'environment(formula)', typically the environment from which 'mv.l1lm' is called.
`subset`	an optional vector specifying a subset of observations to be used in the fitting process.
`na.action`	a function which indicates what should happen when the data contain 'NA's.

Details

The theory behind this function is described in detail in Chapter 13 of the MNM book.

For regular multivariate L2-regression the function lm might be more efficient and offers more methods. Note however that the results given by lm and mv.l1lm may differ slightly due to different divisors of the covariance matrix.

The algorithms for the sign and rank scores are still in an early phase and therefore any feedback is very welcome. For example if p+1 residuals are 0, then the algorithms may not return correct values. Note also that the computations for rank scores might be slow.

Rank regression does not provide an estimate for the intercept parameter is not considered a parameter, a Hodges-Lehmann estimator of the residuals is then an estimate when an interecept term is in the formula. For the one sample case however the function cannot be used for rank scores. We recommend that the regression function should not be used for the one or two sample case. There are distinct functions designed for that purpose. Note furthermore that in the two sample case the covariance matrix returned from the regression function differs slightly from the one returned by the function mv.2sample.est since there matrix A is computed in a different way.

In general it is recommended to use the data argument and specify there the data frame that contains the variables and matrices. For having a matrix Y in a data frame for example the following methods work:

MyData <- data.frame(I(Y),...)

or
MyData <- data.frame(...)
MyData$Y <- Y

Otherwise also the function cbind can be used on the left side of the formula to combine numeric vectors on the fly.

Value

mv.l1ml returns an object of 'class' mvl1lm.

The functions summary is the best choice to view the results. The generic accessor functions coefficients, fitted, residuals and vcov extract various useful features of the value returned by mv.l1ml.

An object of class mv.l1ml is a list wich contains different information depending on the scores and standardization used. To see its content use the function str.

Author(s)

Klaus Nordhausen

References

Oja, H. (2010), Multivariate Nonparametric Methods with R, Springer.

Nordhausen, K. and Oja, H. (2011), Multivariate L1 Methods: The Package MNM, Journal of Statistical Software, 43, 1-28.

Examples

# creating simple data

X <- cbind(rep(1,100),rmvnorm(100,c(0,0,0)) )
B <- matrix(c(4,1,1,0.5,-3,2,2,2),ncol=4, byrow=TRUE)
Y <- X %*% t(B)+ rmvnorm(100,c(0,0), diag(0.2,2))
DAT <- data.frame(x1=X[,2],x2=X[,3], x3=X[,4], Y=I(Y))

# true B
t(B)

# example using identity scores
test1 <- mv.l1lm(Y ~ x1 + x2 + x3, data=DAT)

print(test1)
summary(test1)
coef(test1)
vcov(test1)
head(fitted(test1))
head(residuals(test1))

# example using outer sign scores
test2 <- mv.l1lm(Y ~ x1 + x2 + x3, scores= "s", data=DAT)

print(test2)
summary(test2)
coef(test2)
vcov(test2)
head(fitted(test2))
head(residuals(test2))

# example using inner sign scores
test3 <- mv.l1lm(Y ~ x1 + x2 + x3, scores= "s", stand="i", 
data=DAT)

print(test3)
summary(test3)
coef(test3)
vcov(test3)
head(fitted(test3))
head(residuals(test3))

# example using outer rank scores
test4 <- mv.l1lm(Y ~ x1 + x2 + x3, scores= "r", stand="o", 
data=DAT)

print(test4)
summary(test4)
coef(test4)
vcov(test4)
head(fitted(test4))
head(residuals(test4))

# example using inner rank scores
test5 <- mv.l1lm(Y ~ x1 + x2 + x3, scores= "r", stand="i", 
data=DAT)

print(test5)
summary(test5)
coef(test5)
vcov(test5)
head(fitted(test5))
head(residuals(test5))

# prediction

newData <- data.frame(x1=c(1,-2),x2=c(0.5,0.7), x3=c(-1,-1))
newData
predict(test1,newData)
predict(test2,newData)
predict(test3,newData)
predict(test4,newData)
predict(test5,newData)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

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.

> library(MNM)
Loading required package: ICSNP
Loading required package: mvtnorm
Loading required package: ICS
Loading required package: SpatialNP
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/MNM/mv.l1lm.Rd_%03d_medium.png", width=480, height=480)
> ### Name: mv.l1lm
> ### Title: Linear Regression Based on Identity, Spatial Sign or Spatial
> ###   Rank Scores
> ### Aliases: mv.l1lm
> ### Keywords: multivariate nonparametric regression
> 
> ### ** Examples
> 
> # creating simple data
> 
> X <- cbind(rep(1,100),rmvnorm(100,c(0,0,0)) )
> B <- matrix(c(4,1,1,0.5,-3,2,2,2),ncol=4, byrow=TRUE)
> Y <- X %*% t(B)+ rmvnorm(100,c(0,0), diag(0.2,2))
> DAT <- data.frame(x1=X[,2],x2=X[,3], x3=X[,4], Y=I(Y))
> 
> # true B
> t(B)
     [,1] [,2]
[1,]  4.0   -3
[2,]  1.0    2
[3,]  1.0    2
[4,]  0.5    2
> 
> # example using identity scores
> test1 <- mv.l1lm(Y ~ x1 + x2 + x3, data=DAT)
> 
> print(test1)

Call:
mv.l1lm(formula = Y ~ x1 + x2 + x3, data = DAT)

Coefficients:
             Y1      Y2    
(Intercept)   3.916  -3.039
x1            0.949   1.979
x2            0.955   1.979
x3            0.469   2.046

> summary(test1)

Multivariate regression using identity scores

Call:
mv.l1lm(formula = Y ~ x1 + x2 + x3, data = DAT)

Testing that all coefficients = 0:
Q.2 = 198.1894 with 8 df, p.value < 2.2e-16

Results by response:

Response Y1 :
            Estimate Std. Error
(Intercept)    3.916     0.0436
x1             0.949     0.0398
x2             0.955     0.0419
x3             0.469     0.0432

Response Y2 :
            Estimate Std. Error
(Intercept)    -3.04     0.0398
x1              1.98     0.0364
x2              1.98     0.0383
x3              2.05     0.0395



> coef(test1)
                   Y1        Y2
(Intercept) 3.9159186 -3.038929
x1          0.9490511  1.978639
x2          0.9548634  1.978636
x3          0.4688248  2.046459
> vcov(test1)
               Y1:(Intercept)         Y1:x1         Y1:x2         Y1:x3
Y1:(Intercept)   1.898582e-03  1.088518e-04  1.617530e-04 -4.583422e-05
Y1:x1            1.088518e-04  1.586553e-03 -1.192196e-04 -1.949732e-04
Y1:x2            1.617530e-04 -1.192196e-04  1.753696e-03 -1.326384e-04
Y1:x3           -4.583422e-05 -1.949732e-04 -1.326384e-04  1.863137e-03
Y2:(Intercept)  -1.086264e-04 -6.227896e-06 -9.254612e-06  2.622381e-06
Y2:x1           -6.227896e-06 -9.077378e-05  6.821087e-06  1.115529e-05
Y2:x2           -9.254612e-06  6.821087e-06 -1.003368e-04  7.588836e-06
Y2:x3            2.622381e-06  1.115529e-05  7.588836e-06 -1.065984e-04
               Y2:(Intercept)         Y2:x1         Y2:x2         Y2:x3
Y1:(Intercept)  -1.086264e-04 -6.227896e-06 -9.254612e-06  2.622381e-06
Y1:x1           -6.227896e-06 -9.077378e-05  6.821087e-06  1.115529e-05
Y1:x2           -9.254612e-06  6.821087e-06 -1.003368e-04  7.588836e-06
Y1:x3            2.622381e-06  1.115529e-05  7.588836e-06 -1.065984e-04
Y2:(Intercept)   1.586786e-03  9.097551e-05  1.351890e-04 -3.830707e-05
Y2:x1            9.097551e-05  1.326000e-03 -9.964070e-05 -1.629536e-04
Y2:x2            1.351890e-04 -9.964070e-05  1.465694e-03 -1.108558e-04
Y2:x3           -3.830707e-05 -1.629536e-04 -1.108558e-04  1.557163e-03
> head(fitted(test1))
        Y1        Y2
1 6.473720  4.431287
2 3.571853 -4.640591
3 4.860421 -1.842983
4 6.428832  3.750221
5 2.483841 -6.299658
6 1.832550 -7.222855
> head(residuals(test1))
[1]  0.27742741 -0.08830583  0.01754806 -0.03135157 -0.43719862 -0.36836646
> 
> # example using outer sign scores
> test2 <- mv.l1lm(Y ~ x1 + x2 + x3, scores= "s", data=DAT)
> 
> print(test2)

Call:
mv.l1lm(formula = Y ~ x1 + x2 + x3, scores = "s", data = DAT)

Coefficients:
             Y1      Y2    
(Intercept)   3.931  -3.033
x1            0.946   1.978
x2            0.949   2.011
x3            0.483   2.027

> summary(test2)

Multivariate regression using spatial sign scores and outer standardization

Call:
mv.l1lm(formula = Y ~ x1 + x2 + x3, scores = "s", data = DAT)

Testing that all coefficients = 0:
Q.2 = 190.2277 with 8 df, p.value < 2.2e-16

Results by response:

Response Y1 :
            Estimate Std. Error
(Intercept)    3.931    0.00148
x1             0.946    0.00135
x2             0.949    0.00142
x3             0.483    0.00146

Response Y2 :
            Estimate Std. Error
(Intercept)    -3.03     0.0381
x1              1.98     0.0348
x2              2.01     0.0366
x3              2.03     0.0378



> coef(test2)
                   Y1        Y2
(Intercept) 3.9306725 -3.032808
x1          0.9458760  1.978282
x2          0.9485644  2.011111
x3          0.4831106  2.027273
> vcov(test2)
               Y1:(Intercept)         Y1:x1         Y1:x2         Y1:x3
Y1:(Intercept)   2.183849e-06  1.252070e-07  1.860568e-07 -5.272094e-08
Y1:x1            1.252070e-07  1.824937e-06 -1.371327e-07 -2.242684e-07
Y1:x2            1.860568e-07 -1.371327e-07  2.017193e-06 -1.525677e-07
Y1:x3           -5.272094e-08 -2.242684e-07 -1.525677e-07  2.143079e-06
Y2:(Intercept)   5.590279e-05  3.205085e-06  4.762735e-06 -1.349566e-06
Y2:x1            3.205085e-06  4.671524e-05 -3.510361e-06 -5.740887e-06
Y2:x2            4.762735e-06 -3.510361e-06  5.163668e-05 -3.905471e-06
Y2:x3           -1.349566e-06 -5.740887e-06 -3.905471e-06  5.485914e-05
               Y2:(Intercept)         Y2:x1         Y2:x2         Y2:x3
Y1:(Intercept)   5.590279e-05  3.205085e-06  4.762735e-06 -1.349566e-06
Y1:x1            3.205085e-06  4.671524e-05 -3.510361e-06 -5.740887e-06
Y1:x2            4.762735e-06 -3.510361e-06  5.163668e-05 -3.905471e-06
Y1:x3           -1.349566e-06 -5.740887e-06 -3.905471e-06  5.485914e-05
Y2:(Intercept)   1.453221e-03  8.331777e-05  1.238096e-04 -3.508262e-05
Y2:x1            8.331777e-05  1.214386e-03 -9.125357e-05 -1.492372e-04
Y2:x2            1.238096e-04 -9.125357e-05  1.342321e-03 -1.015246e-04
Y2:x3           -3.508262e-05 -1.492372e-04 -1.015246e-04  1.426090e-03
> head(fitted(test2))
        Y1        Y2
1 6.513339  4.379235
2 3.573541 -4.606326
3 4.856783 -1.781474
4 6.455246  3.760420
5 2.504209 -6.342775
6 1.861075 -7.268233
> head(residuals(test2))
[1]  0.23780796 -0.08999386  0.02118659 -0.05776513 -0.45756641 -0.39689101
> 
> # example using inner sign scores
> test3 <- mv.l1lm(Y ~ x1 + x2 + x3, scores= "s", stand="i", 
+ data=DAT)
> 
> print(test3)

Call:
mv.l1lm(formula = Y ~ x1 + x2 + x3, scores = "s", stand = "i",     data = DAT)

Coefficients:
             Y1      Y2    
(Intercept)   3.934  -3.033
x1            0.948   1.977
x2            0.948   2.010
x3            0.487   2.027

> summary(test3)

Multivariate regression using spatial sign scores and inner standardization

Call:
mv.l1lm(formula = Y ~ x1 + x2 + x3, scores = "s", stand = "i",     data = DAT)

Testing that all coefficients = 0:
Q.2 = 178.9022 with 8 df, p.value < 2.2e-16

Results by response:

Response Y1 :
            Estimate Std. Error
(Intercept)    3.934    0.00409
x1             0.948    0.00374
x2             0.948    0.00394
x3             0.487    0.00406

Response Y2 :
            Estimate Std. Error
(Intercept)    -3.03     0.0363
x1              1.98     0.0332
x2              2.01     0.0349
x3              2.03     0.0359



> coef(test3)
                   Y1        Y2
(Intercept) 3.9343050 -3.033205
x1          0.9479176  1.977361
x2          0.9479922  2.010450
x3          0.4869700  2.027055
> vcov(test3)
              :(Intercept)           :x1           :x2           :x3
:(Intercept)  1.676565e-05  9.612279e-07  1.428378e-06 -4.047443e-07
:x1           9.612279e-07  1.401024e-05 -1.052782e-06 -1.721733e-06
:x2           1.428378e-06 -1.052782e-06  1.548621e-05 -1.171279e-06
:x3          -4.047443e-07 -1.721733e-06 -1.171279e-06  1.645265e-05
:(Intercept)  1.485238e-04  8.515341e-06  1.265374e-05 -3.585555e-06
:x1           8.515341e-06  1.241141e-04 -9.326405e-06 -1.525252e-05
:x2           1.265374e-05 -9.326405e-06  1.371895e-04 -1.037614e-05
:x3          -3.585555e-06 -1.525252e-05 -1.037614e-05  1.457510e-04
              :(Intercept)           :x1           :x2           :x3
:(Intercept)  1.485238e-04  8.515341e-06  1.265374e-05 -3.585555e-06
:x1           8.515341e-06  1.241141e-04 -9.326405e-06 -1.525252e-05
:x2           1.265374e-05 -9.326405e-06  1.371895e-04 -1.037614e-05
:x3          -3.585555e-06 -1.525252e-05 -1.037614e-05  1.457510e-04
:(Intercept)  1.316660e-03  7.548833e-05  1.121752e-04 -3.178588e-05
:x1           7.548833e-05  1.100269e-03 -8.267840e-05 -1.352133e-04
:x2           1.121752e-04 -8.267840e-05  1.216182e-03 -9.198429e-05
:x3          -3.178588e-05 -1.352133e-04 -9.198429e-05  1.292080e-03
> head(fitted(test3))
        Y1        Y2
1 6.529702  4.376672
2 3.573098 -4.606483
3 4.857045 -1.782613
4 6.465826  3.758195
5 2.508374 -6.342291
6 1.864583 -7.266978
> head(residuals(test3))
[1]  0.22144556 -0.08955058  0.02092388 -0.06834594 -0.46173220 -0.40039947
> 
> # example using outer rank scores
> test4 <- mv.l1lm(Y ~ x1 + x2 + x3, scores= "r", stand="o", 
+ data=DAT)
> 
> print(test4)

Call:
mv.l1lm(formula = Y ~ x1 + x2 + x3, scores = "r", stand = "o",     data = DAT)

Coefficients:
    Y1     Y2   
x1  0.946  1.975
x2  0.956  1.996
x3  0.480  2.040

> summary(test4)

Multivariate regression using spatial rank scores and outer standardization

Call:
mv.l1lm(formula = Y ~ x1 + x2 + x3, scores = "r", stand = "o",     data = DAT)

Outer HL-estimator for the residuals (intercept):
              Y1    Y2
(Intercept) 3.92 -3.04

Testing that all coefficients = 0:
Q.2 = 140.0406 with 6 df, p.value < 2.2e-16

Results by response:

Response Y1 :
   Estimate Std. Error
x1    0.946     0.0397
x2    0.956     0.0417
x3    0.480     0.0430

Response Y2 :
   Estimate Std. Error
x1     1.98     0.0359
x2     2.00     0.0378
x3     2.04     0.0390



> coef(test4)
          Y1       Y2
x1 0.9458646 1.975473
x2 0.9557754 1.996265
x3 0.4804507 2.039694
> vcov(test4)
              Y1:x1         Y1:x2         Y1:x3         Y2:x1         Y2:x2
Y1:x1  0.0015739642 -1.182737e-04 -1.934262e-04 -0.0001667949  1.253360e-05
Y1:x2 -0.0001182737  1.739781e-03 -1.315860e-04  0.0000125336 -1.843667e-04
Y1:x3 -0.0001934262 -1.315860e-04  1.848354e-03  0.0000204976  1.394432e-05
Y2:x1 -0.0001667949  1.253360e-05  2.049760e-05  0.0012923221 -9.711000e-05
Y2:x2  0.0000125336 -1.843667e-04  1.394432e-05 -0.0000971100  1.428468e-03
Y2:x3  0.0000204976  1.394432e-05 -1.958723e-04 -0.0001588149 -1.080402e-04
              Y2:x3
Y1:x1  2.049760e-05
Y1:x2  1.394432e-05
Y1:x3 -1.958723e-04
Y2:x1 -1.588149e-04
Y2:x2 -1.080402e-04
Y2:x3  1.517613e-03
> head(fitted(test4))
        Y1        Y2
1 6.497539  4.397881
2 3.572222 -4.629577
3 4.861722 -1.817906
4 6.452323  3.752776
5 2.486625 -6.330584
6 1.843572 -7.250399
> head(residuals(test4))
[1]  0.25360876 -0.08867427  0.01624743 -0.05484233 -0.43998319 -0.37938821
> 
> # example using inner rank scores
> test5 <- mv.l1lm(Y ~ x1 + x2 + x3, scores= "r", stand="i", 
+ data=DAT)
> 
> print(test5)

Call:
mv.l1lm(formula = Y ~ x1 + x2 + x3, scores = "r", stand = "i",     data = DAT)

Coefficients:
    Y1     Y2   
x1  0.946  1.975
x2  0.956  1.996
x3  0.481  2.038

> summary(test5)

Multivariate regression using spatial rank scores and inner standardization

Call:
mv.l1lm(formula = Y ~ x1 + x2 + x3, scores = "r", stand = "i",     data = DAT)

Inner HL-estimator for the residuals (intercept):
              Y1    Y2
(Intercept) 3.92 -3.04

Testing that all coefficients = 0:
Q.2 = 139.9164 with 6 df, p.value < 2.2e-16

Results by response:

Response Y1 :
   Estimate Std. Error
x1    0.946     0.0397
x2    0.956     0.0417
x3    0.481     0.0430

Response Y2 :
   Estimate Std. Error
x1     1.97     0.0364
x2     2.00     0.0382
x3     2.04     0.0394



> coef(test5)
          Y1       Y2
x1 0.9462485 1.974523
x2 0.9555783 1.996343
x3 0.4814340 2.038414
> vcov(test5)
              Y1:x1         Y1:x2         Y1:x3         Y2:x1         Y2:x2
Y1:x1  1.574548e-03 -1.183175e-04 -1.934979e-04 -1.797794e-04  1.350931e-05
Y1:x2 -1.183175e-04  1.740426e-03 -1.316348e-04  1.350931e-05 -1.987191e-04
Y1:x3 -1.934979e-04 -1.316348e-04  1.849040e-03  2.209328e-05  1.502985e-05
Y2:x1 -1.797794e-04  1.350931e-05  2.209328e-05  1.322292e-03 -9.936207e-05
Y2:x2  1.350931e-05 -1.987191e-04  1.502985e-05 -9.936207e-05  1.461595e-03
Y2:x3  2.209328e-05  1.502985e-05 -2.111204e-04 -1.624979e-04 -1.105458e-04
              Y2:x3
Y1:x1  2.209328e-05
Y1:x2  1.502985e-05
Y1:x3 -2.111204e-04
Y2:x1 -1.624979e-04
Y2:x2 -1.105458e-04
Y2:x3  1.552808e-03
> head(fitted(test5))
        Y1        Y2
1 6.501033  4.392584
2 3.571745 -4.628689
3 4.861325 -1.817460
4 6.454451  3.749595
5 2.487340 -6.331172
6 1.844253 -7.250500
> head(residuals(test5))
[1]  0.25011420 -0.08819719  0.01664389 -0.05697079 -0.44069731 -0.38006966
> 
> # prediction
> 
> newData <- data.frame(x1=c(1,-2),x2=c(0.5,0.7), x3=c(-1,-1))
> newData
  x1  x2 x3
1  1 0.5 -1
2 -2 0.7 -1
> predict(test1,newData)
        Y1        Y2
1 4.873577 -2.117431
2 2.217396 -7.657622
> predict(test2,newData)
        Y1        Y2
1 4.867720 -2.076244
2 2.219805 -7.608868
> predict(test3,newData)
        Y1        Y2
1 4.869249 -2.077674
2 2.215094 -7.607667
> predict(test4,newData)
        Y1        Y2
1 4.867745 -2.107274
2 2.221306 -7.634439
> predict(test5,newData)
        Y1        Y2
1 4.867582 -2.107418
2 2.219952 -7.631718
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>