R/lsfitNci.R

lsfitNci <-
function(x,y,alpha=.05){
#
# Compute confidence for least squares
# regression using heteroscedastic method
# recommended by Long and Ervin (2000).
#
x<-as.matrix(x)
if(nrow(x) != length(y))stop("Length of y does not match number of x values")
m<-cbind(x,y)
m<-elimna(m)
y<-m[,ncol(x)+1]
temp<-lsfit(x,y)
x<-cbind(rep(1,nrow(x)),m[,1:ncol(x)])
xtx<-solve(t(x)%*%x)
h<-diag(x%*%xtx%*%t(x))
hc3<-xtx%*%t(x)%*%diag(temp$res^2/(1-h)^2)%*%x%*%xtx
df<-nrow(x)-ncol(x)
crit<-qt(1-alpha/2,df)
al<-ncol(x)
ci<-matrix(NA,nrow=al,ncol=3)
for(j in 1:al){
ci[j,1]<-j
ci[j,2]<-temp$coef[j]-crit*sqrt(hc3[j,j])
ci[j,3]<-temp$coef[j]+crit*sqrt(hc3[j,j])
}
print("Confidence intervals for intercept followed by slopes:")
list(ci=ci,stand.errors=sqrt(diag(hc3)))
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.