https://github.com/cran/robCompositions
Tip revision: bea0f6203c78ee704dca8deb7c74cc48e13ee950 authored by Matthias Templ on 15 May 2012, 15:35:41 UTC
version 1.6.0
version 1.6.0
Tip revision: bea0f62
lmCoDaX.R
lmCoDaX <- function(y, X, method = "robust"){
ilrregression <- function(X,y){
# delivers appropriate inference for regression of y on a compositional matrix X
# PF, Aug 18, 2010
#
# classical regression
d <- data.frame(y=y,X=X)
lmcla <- lm(y~.,data=d)
lmcla.sum <- summary(lmcla)
# ilr regressions:
require(robCompositions)
ilr.sum <- lmcla.sum
for (j in 1:ncol(X)){
Zj <- -robCompositions::ilr(cbind(X[,j],X[,-j]))
dj <- data.frame(y=y,Z=Zj)
res <- lm(y~.,data=dj)
res.sum <- summary(res)
if (j==1){
ilr.sum$coefficients[1:2,] <- res.sum$coefficients[1:2,]
ilr.sum$residuals <- res.sum$residuals
ilr.sum$sigma <- res.sum$sigma
ilr.sum$r.squared <- res.sum$r.squared
ilr.sum$adj.r.squared <- res.sum$adj.r.squared
ilr.sum$fstatistic <- res.sum$fstatistic
}
else{
ilr.sum$coefficients[j+1,] <- res.sum$coefficients[2,]
}
}
list(lm = lmcla, lm=lmcla.sum,ilr=ilr.sum)
}
robilrregression <- function(X,y){
# delivers appropriate inference for robust regression of y on a compositional matrix X
# PF, Aug 18, 2010
#
# classical regression
require(robustbase)
d <- data.frame(y=y,X=X)
lmcla <- ltsReg(y~.,data=d)
lmcla.sum <- summary(lmcla)
# ilr regressions:
require(robCompositions)
ilr.sum <- lmcla.sum
for (j in 1:ncol(X)){
Zj <- -robCompositions::ilr(cbind(X[,j],X[,-j]))
dj <- data.frame(y=y,Z=Zj)
res <- ltsReg(y~.,data=dj)
res.sum <- summary(res)
if (j==1){
ilr.sum$coefficients[1:2,] <- res.sum$coefficients[1:2,]
ilr.sum$residuals <- res.sum$residuals
ilr.sum$sigma <- res.sum$sigma
ilr.sum$r.squared <- res.sum$r.squared
ilr.sum$adj.r.squared <- res.sum$adj.r.squared
ilr.sum$fstatistic <- res.sum$fstatistic
}
else{
ilr.sum$coefficients[j+1,] <- res.sum$coefficients[2,]
}
}
list(lm = lmcla, lm=lmcla.sum,ilr=ilr.sum)
}
if( method =="classical"){
reg <- ilrregression(X,y)
} else if(method=="robust"){
reg <- robilrregression(X,y)
}
return(reg)
}