If you try to fit a regression model with p>n parameter - what happens?
In [1]:
library(MASS)
p<-200
n<-100
k<-5
XX<-mvrnorm(n,rep(0,p),diag(rep(1,p)))
y<-XX[,1:k]%*%rep(1,k)+rnorm(n)
mm<-lm(y~XX-1) #without the intercept
summary(mm)
plot(mm$coef,xlab="coefficient",ylab="estimate")
abline(v=k)
Call:
lm(formula = y ~ XX - 1)
Residuals:
ALL 100 residuals are 0: no residual degrees of freedom!
Coefficients: (100 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
XX1 0.89676 NaN NaN NaN
XX2 0.93325 NaN NaN NaN
XX3 2.07332 NaN NaN NaN
XX4 2.41896 NaN NaN NaN
XX5 -3.71722 NaN NaN NaN
XX6 0.41146 NaN NaN NaN
XX7 -0.12414 NaN NaN NaN
XX8 1.86365 NaN NaN NaN
XX9 -0.93305 NaN NaN NaN
XX10 1.44008 NaN NaN NaN
XX11 -3.32781 NaN NaN NaN
XX12 1.99823 NaN NaN NaN
XX13 -1.02657 NaN NaN NaN
XX14 -0.24226 NaN NaN NaN
XX15 2.35035 NaN NaN NaN
XX16 -0.60457 NaN NaN NaN
XX17 -2.02729 NaN NaN NaN
XX18 0.04420 NaN NaN NaN
XX19 2.92674 NaN NaN NaN
XX20 1.71296 NaN NaN NaN
XX21 -2.72378 NaN NaN NaN
XX22 0.52028 NaN NaN NaN
XX23 1.01874 NaN NaN NaN
XX24 -0.17535 NaN NaN NaN
XX25 -3.28576 NaN NaN NaN
XX26 1.18211 NaN NaN NaN
XX27 -0.89578 NaN NaN NaN
XX28 -0.49626 NaN NaN NaN
XX29 0.14176 NaN NaN NaN
XX30 5.01389 NaN NaN NaN
XX31 0.71798 NaN NaN NaN
XX32 0.46278 NaN NaN NaN
XX33 2.34123 NaN NaN NaN
XX34 3.98623 NaN NaN NaN
XX35 -2.77560 NaN NaN NaN
XX36 -2.00876 NaN NaN NaN
XX37 1.49524 NaN NaN NaN
XX38 -0.02343 NaN NaN NaN
XX39 -0.66684 NaN NaN NaN
XX40 3.04627 NaN NaN NaN
XX41 -2.75145 NaN NaN NaN
XX42 -4.11748 NaN NaN NaN
XX43 -3.71846 NaN NaN NaN
XX44 0.16381 NaN NaN NaN
XX45 5.50446 NaN NaN NaN
XX46 -0.54734 NaN NaN NaN
XX47 4.81076 NaN NaN NaN
XX48 2.35173 NaN NaN NaN
XX49 1.61631 NaN NaN NaN
XX50 -1.20593 NaN NaN NaN
XX51 1.29742 NaN NaN NaN
XX52 2.86436 NaN NaN NaN
XX53 0.97924 NaN NaN NaN
XX54 -0.34676 NaN NaN NaN
XX55 -0.77451 NaN NaN NaN
XX56 3.35063 NaN NaN NaN
XX57 -1.50725 NaN NaN NaN
XX58 -3.10193 NaN NaN NaN
XX59 2.26055 NaN NaN NaN
XX60 -1.65352 NaN NaN NaN
XX61 4.05446 NaN NaN NaN
XX62 -0.09644 NaN NaN NaN
XX63 1.10096 NaN NaN NaN
XX64 1.00321 NaN NaN NaN
XX65 1.51325 NaN NaN NaN
XX66 5.34817 NaN NaN NaN
XX67 -4.03162 NaN NaN NaN
XX68 2.98884 NaN NaN NaN
XX69 -3.04321 NaN NaN NaN
XX70 0.54842 NaN NaN NaN
XX71 -1.02863 NaN NaN NaN
XX72 -0.55798 NaN NaN NaN
XX73 2.79298 NaN NaN NaN
XX74 2.46611 NaN NaN NaN
XX75 -0.37112 NaN NaN NaN
XX76 -2.96653 NaN NaN NaN
XX77 2.91184 NaN NaN NaN
XX78 -2.31267 NaN NaN NaN
XX79 0.04920 NaN NaN NaN
XX80 4.90122 NaN NaN NaN
XX81 -1.37344 NaN NaN NaN
XX82 0.77549 NaN NaN NaN
XX83 1.25408 NaN NaN NaN
XX84 -2.28155 NaN NaN NaN
XX85 -1.15282 NaN NaN NaN
XX86 -2.82518 NaN NaN NaN
XX87 0.35251 NaN NaN NaN
XX88 0.30585 NaN NaN NaN
XX89 3.09787 NaN NaN NaN
XX90 2.27595 NaN NaN NaN
XX91 -1.36101 NaN NaN NaN
XX92 6.85997 NaN NaN NaN
XX93 0.49749 NaN NaN NaN
XX94 -2.35413 NaN NaN NaN
XX95 2.22426 NaN NaN NaN
XX96 -2.52583 NaN NaN NaN
XX97 2.04452 NaN NaN NaN
XX98 -1.25311 NaN NaN NaN
XX99 -2.85343 NaN NaN NaN
XX100 1.65932 NaN NaN NaN
XX101 NA NA NA NA
XX102 NA NA NA NA
XX103 NA NA NA NA
XX104 NA NA NA NA
XX105 NA NA NA NA
XX106 NA NA NA NA
XX107 NA NA NA NA
XX108 NA NA NA NA
XX109 NA NA NA NA
XX110 NA NA NA NA
XX111 NA NA NA NA
XX112 NA NA NA NA
XX113 NA NA NA NA
XX114 NA NA NA NA
XX115 NA NA NA NA
XX116 NA NA NA NA
XX117 NA NA NA NA
XX118 NA NA NA NA
XX119 NA NA NA NA
XX120 NA NA NA NA
XX121 NA NA NA NA
XX122 NA NA NA NA
XX123 NA NA NA NA
XX124 NA NA NA NA
XX125 NA NA NA NA
XX126 NA NA NA NA
XX127 NA NA NA NA
XX128 NA NA NA NA
XX129 NA NA NA NA
XX130 NA NA NA NA
XX131 NA NA NA NA
XX132 NA NA NA NA
XX133 NA NA NA NA
XX134 NA NA NA NA
XX135 NA NA NA NA
XX136 NA NA NA NA
XX137 NA NA NA NA
XX138 NA NA NA NA
XX139 NA NA NA NA
XX140 NA NA NA NA
XX141 NA NA NA NA
XX142 NA NA NA NA
XX143 NA NA NA NA
XX144 NA NA NA NA
XX145 NA NA NA NA
XX146 NA NA NA NA
XX147 NA NA NA NA
XX148 NA NA NA NA
XX149 NA NA NA NA
XX150 NA NA NA NA
XX151 NA NA NA NA
XX152 NA NA NA NA
XX153 NA NA NA NA
XX154 NA NA NA NA
XX155 NA NA NA NA
XX156 NA NA NA NA
XX157 NA NA NA NA
XX158 NA NA NA NA
XX159 NA NA NA NA
XX160 NA NA NA NA
XX161 NA NA NA NA
XX162 NA NA NA NA
XX163 NA NA NA NA
XX164 NA NA NA NA
XX165 NA NA NA NA
XX166 NA NA NA NA
XX167 NA NA NA NA
XX168 NA NA NA NA
XX169 NA NA NA NA
XX170 NA NA NA NA
XX171 NA NA NA NA
XX172 NA NA NA NA
XX173 NA NA NA NA
XX174 NA NA NA NA
XX175 NA NA NA NA
XX176 NA NA NA NA
XX177 NA NA NA NA
XX178 NA NA NA NA
XX179 NA NA NA NA
XX180 NA NA NA NA
XX181 NA NA NA NA
XX182 NA NA NA NA
XX183 NA NA NA NA
XX184 NA NA NA NA
XX185 NA NA NA NA
XX186 NA NA NA NA
XX187 NA NA NA NA
XX188 NA NA NA NA
XX189 NA NA NA NA
XX190 NA NA NA NA
XX191 NA NA NA NA
XX192 NA NA NA NA
XX193 NA NA NA NA
XX194 NA NA NA NA
XX195 NA NA NA NA
XX196 NA NA NA NA
XX197 NA NA NA NA
XX198 NA NA NA NA
XX199 NA NA NA NA
XX200 NA NA NA NA
Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: NaN
F-statistic: NaN on 100 and 0 DF, p-value: NA
Did NOT work. The first n coefficients are estimated (quite arbitrary!) and the rest are not. Moreover, you used up all your degrees of freedom by doing estimation so do not obtain any estimates for the noise variance, $\sigma$.
Try different values of p at home. Even if you could fit the model, are the estimates any good when p is close to n?
What if we use ridge-regression?
In [2]:
p<-200
n<-100
k<-5
lambda<-.2
XX<-mvrnorm(n,rep(0,p),diag(rep(1,p)))
y<-XX[,1:k]%*%rep(1,k)+rnorm(n)
bb<-solve(t(XX)%*%XX+lambda*diag(rep(1,p)))%*%(t(XX)%*%y) #without the intercept
plot(bb,xlab="coefficient",ylab="estimate")
abline(v=k)
This regularization helped reduce the numerical instability of the inverse of X'X.
Can we achieve regularization through other means? Let's try using gradient descent - starting with all coefficients equal to 0.
In [12]:
p<-200
n<-100
k<-5
Niter<-500 # number of iterations
eta<-.05 # learning rate
lambda<-0 # using ridge penalty or not
BETA<-matrix(0,200,Niter)
for (zz in (2:Niter)) {
BETA[,zz] <- BETA[,zz-1] + eta*(t(XX)%*%(y-XX%*%BETA[,zz-1])+lambda*BETA[,zz-1])/n
}
plot(seq(1,Niter),BETA[1,],col=1,ylim=c(-.25,1.25))
for (z2 in (2:p)) {
lines(seq(1,Niter),BETA[z2,],col=z2,lty=z2)
}
plot(BETA[,Niter])
abline(v=k)
