机器算法验证 - 基于簇夹心 VCV 的惩罚样条置信区间 - 吾爱随笔录

这是我在这里的第一篇文章，但我从这个论坛在谷歌搜索结果中弹出的结果中受益匪浅。

我一直在使用惩罚样条自学半参数回归。基本上，

\hat{y} = X (X^{'} X + λ D)^{- 1} X^{'} y

$\hat{y} = X(X'X + \lambda D)^{-1}X'y$ 除了我的原始变量之外，X 还充满了结，

λ

$\lambda$ 是应用于结的惩罚项。（如果您想了解这些内容，请参阅 Wood 2006 或 Ruppert, Wand & Carrol 2003）。

我的数据具有聚类自相关，并且可能也具有异方差性。我的培训是计量经济学，所以我自然想将方差-协方差矩阵表示为

X (X^{'} X + λ D)^{- 1} (\sum_{j} X_{j}^{'} ϵ_{j} ϵ_{j}^{'} X) (X^{'} X + λ D)^{- 1} X^{'}

$X(X'X + \lambda D)^{-1}\left(\displaystyle\sum_j X'_j \epsilon_j \epsilon_j'X\right)(X'X + \lambda D)^{-1}X'$ 以上假设集群之间是独立的，但不是在它们内部，并且不假设同方差。这应该给我正确的se。

我的问题是我不只关心se。使用半参数，我需要图上的置信区间。

一段时间后编辑
：我已经开发了一种通过 hacking 来做到这一点的准蛮力方式mgcv。一般的做法是计算我自己的参数协方差矩阵并将其填充到mgcv's. 然后其他函数使用它来计算 CI 并在图上给出 CI 这是一项正在进行中的工作，我对它还没有信心，但我把它提出来，希望它可能是迈向某事的第一步这很有效，可以邀请反馈，并可能最终对其他人有用。这是目前的情况：

另一个编辑：该函数现在尝试通过反转参数协方差矩阵和减法来获得完整的惩罚矩阵 $X'X$ . 奇怪的是，由提供的参数协方差矩阵mgcv通常是不可逆的。在删除自动丢弃的冗余假人之后 ，我不知道为什么会出现这种情况。有人在这里有什么想法吗？ 我目前正在使用广义逆作为一种技巧，我不确定究竟会增加什么样的失真。

clusterGAM = function(m,clustvar){
    require(MASS)
    mm = predict(m,type='lpmatrix') #Get model X matrix
    mm = mm[,colnames(mm) %in% names(coefficients(m))[rowSums(m$Vp)!=0]] #Cut out coefficients from model matrix that get auto-dropped due to collinearity
	meat = matrix(rep(0,ncol(mm)*ncol(mm)),ncol(mm)) #Pre-fill a meat matrix
	for (i in 1:length(unique(clustvar))){#This sums over the clusters
		X = mm[clustvar == unique(clustvar)[i],]
		e = residuals(m)[clustvar == unique(clustvar)[i]]
		meat = meat+t(X) %*% e %*% t(e) %*% X
		}
	M <- length(unique(clustvar))
	N <- length(clustvar)
	K <- sum(m$edf) #use EDF instead of K
    dfc <- (M/(M-1))*((N-1)/(N-K))#Degrees of freedom correction
    Vp = vcov(m,dispersion=1)
    Vp = Vp[rowSums(m$Vp)!=0,rowSums(m$Vp)!=0]

    S = tryCatch(solve(Vp) - t(mm)%*%mm, error=function(e) e, finally=NA)
    if(inherits(S, "error")) {
        S = ginv(Vp) - t(mm)%*%mm
        print("Warning:  Generalized inverse used to extract the penalty matrix")
        }   
    vcv = tryCatch(dfc*solve(t(mm)%*%mm+S)%*%meat %*% solve(t(mm)%*%mm+S), error=function(e) e, finally=NA)
    if(inherits(vcv, "error")) {
        vcv = dfc*ginv(t(mm)%*%mm+S)%*%meat %*% ginv(t(mm)%*%mm+S)
        print("Warning:  Generalized inverse used to calculate the Vp matrix")
        } #Get the vcv
    if (length(which(rowSums(m$Vp)==0)!=0)){	#If any zero rows got knocked out, put them back in
		zeros = which(rowSums(m$Vp)==0)
		for (i in 1:length(zeros)){
			vcv = rbind(vcv[0:(zeros[i]-1),],0,vcv[(zeros[i]):nrow(vcv),])
			vcv = cbind(vcv[,0:(zeros[i]-1)],0,vcv[,(zeros[i]):ncol(vcv)])
			}
		}
	m$Vp = vcv #stuff the new vcv into the place where the plots and the summary gets made from
    m
    }

*编辑为能够处理张量

一些数据来运行它

load(broken link)
head(sdat)
N = nrow(sdat)
#set.seed(417)
sdat$x = with(sdat,log(y+1+runif(N)*T)+rexp(N)) #Making a fake smooth covariate
    plot(sdat$x,sdat$y)
    m = gam(y~ID+G:t+T+s(x)-1,data=sdat)#LSDV panel model with time x group FEs and one smooth function
    par(mfrow=c(1,2))
    plot(m,ylim=c(-.75,2))
    summary(m)
    m = clusterGAM(m,sdat$G)
plot(m,ylim=c(-.75,2)) #the confidence intervals shrink a bit.  This seems counterintuitive.  
summary(m) #SE on T gets a little bit bigger, as it probably should if you account for repeated observations

#Now try it with a tensor
sdat$z = with(sdat,y+x+rnorm(N))
m = gam(y~ID+G:t+T+te(x,z)-1,data=sdat)
summary(m)
vis.gam(m,view=c('x','z'),se=2,theta=150,phi=20)
m = clusterGAM(m,sdat$G)
vis.gam(m,view=c('x','z'),se=2,theta=150,phi=20)

它可能被我对所涉及的数学的不完全理解或mgcv. 就目前而言，参数假人的标准误差以看似合理的方式膨胀，但平滑项的置信区间似乎缩小了，这似乎是错误的。

任何关于改进这一点的意见将不胜感激。

最后一点：我不是统计学博士。我是一名从事应用工作的研究生，可能花费的时间比我应该让自己被功能形式假设所困扰的时间多得多。