0
$\begingroup$
n = 50
set.seed(100)
x = matrix(runif(n, -2, 2), nrow=n)
y = 2 + 0.75*sin(x) - 0.75*cos(x) + rnorm(n, 0, 0.2)


nknots = 12
mat.tpb = function(x, nknots){
    n = length(x)
    knots = seq(-2, 2, len=nknots)
    zb = cbind(1, x ,x^2, x^3)
    zt = matrix(NA, nrow=n, ncol=nknots)
    for(i in 1:nknots){zt[,i] = ifelse(x>knots[i], (x-knots[i]) ^3, 0)}
    cbind(zb, zt)}

z = mat.tpb(x, nknots)

I wanted to do this:

> summary(lm(y~z))

Call:
lm(formula = y ~ z)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.35945 -0.06861  0.00875  0.08329  0.39890 

Coefficients: (3 not defined because of singularities)
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) 30636.8811 12890.7059   2.377  0.02291 * 
z1                  NA         NA      NA       NA   
z2          55942.6648 23543.5397   2.376  0.02294 * 
z3          34019.1008 14319.8448   2.376  0.02296 * 
z4           6887.7980  2899.9304   2.375  0.02299 * 
z5                  NA         NA      NA       NA   
z6          -7099.3431  2983.7929  -2.379  0.02277 * 
z7            243.7723    94.1488   2.589  0.01379 * 
z8            -46.0888    15.2702  -3.018  0.00465 **
z9             22.0511     9.1779   2.403  0.02156 * 
z10           -16.0403     8.5338  -1.880  0.06827 . 
z11            13.6275     8.5711   1.590  0.12059   
z12            -6.5516     9.1134  -0.719  0.47685   
z13            -3.5033    10.1582  -0.345  0.73220   
z14            10.1460    12.0295   0.843  0.40456   
z15             0.9644    27.5610   0.035  0.97228   
z16                 NA         NA      NA       NA   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1923 on 36 degrees of freedom
Multiple R-squared:  0.951, Adjusted R-squared:  0.9332 
F-statistic: 53.69 on 13 and 36 DF,  p-value: < 2.2e-16

But it seems the design matrix z is singular,

so I used Moore-Penrose generalized inverse matrix to do this:

> MASS::ginv(t(z)%*%z)%*%t(z)%*%y
             [,1]
 [1,]  0.43424264
 [2,] -0.45336692
 [3,]  0.28938730
 [4,]  0.17283360
 [5,] -0.05730462
 [6,] -0.11531967
 [7,]  1.12054017
 [8,]  0.30900705
 [9,] -3.79858631
[10,]  1.52160075
[11,]  3.67809815
[12,] -4.15539087
[13,] -0.64600815
[14,]  5.77855917
[15,]  1.53625681
[16,]  0.00000000

But I am not sure whether the calculation is right, which means I do not know the calculation of regression coefficient calculation when using the Moore-Penrose generalized inverse matrix.

Please tell me whether the regression coefficient calculation by Moore-Penrose generalized inverse matrix is correct or not.

Thank you very much.

$\endgroup$
  • $\begingroup$ If the matrix is full rank, the Moore-Penrose inverse equals the "regular" inverse. So in that sense using Moore-Penrose is correct. Note, however, that there are alternatives for the case $p > n$ or singular design matrices. In the latter case, figuring out why your design matrix is singular is likely to be far more productive than using a generalized inverse; alternatively, sparsity-inducing procedures such as Lasso can be used. $\endgroup$ – jbowman Dec 1 '19 at 15:56
  • $\begingroup$ It's obviously incorrect because it includes only 16 variables rather than the 17 you used for lm. The very first thing you need to do is reduce your example to something tiny and manageable--consider, say, a dataset with one variable and two or three observations. It's unfair to ask your readers to wade through the code and the output for such a large problem when, with a little work on your part, something much simpler and accessible will do. Moreover, it's highly likely that in the process of creating such a minimal reproducible example you will discover the answer yourself. $\endgroup$ – whuber Dec 1 '19 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.