### R code from vignette source 'Chapter05_LinearModel.prnw' ################################################### ### code chunk number 1: Chapter05_LinearModel.prnw:34-35 ################################################### options(show.signif.stars=FALSE) ################################################### ### code chunk number 2: Chapter05_LinearModel.prnw:136-140 ################################################### data(Orange) names(Orange) model = lm(circumference~age,data=Orange) model ################################################### ### code chunk number 3: Chapter05_LinearModel.prnw:150-153 ################################################### plot(circumference~age,data=Orange) abline(model) points(predict(model)~age,data=Orange,pch=20) ################################################### ### code chunk number 4: Chapter05_LinearModel.prnw:157-159 ################################################### var(model) summary(model)$sigma^2 ################################################### ### code chunk number 5: Chapter05_LinearModel.prnw:163-170 ################################################### #confidence <- function(model,conf=0.95,alpha=1-conf) { #% tvalue <- qt(c(lower=alpha/2,upper=1-alpha/2),df=model$df.residual) # Delta <- outer(sqrt(diag(vcov(model))) ,tvalue) # coef(model)+Delta #} #confidence(model) confint(model) ################################################### ### code chunk number 6: Chapter05_LinearModel.prnw:179-185 ################################################### mynew = seq(from=min(Orange$age),to=max(Orange$age),length.out=100) mynew = data.frame(age=mynew) options(max.print=40) (CI = predict(model,newdata=mynew, interval="confidence")) (PI = predict(model,newdata=mynew, interval="prediction")) options(max.print=99999) ################################################### ### code chunk number 7: Chapter05_LinearModel.prnw:194-196 ################################################### matlines(mynew, CI, col="red", lty=2) matlines(mynew, PI, col="black", lty=3) ################################################### ### code chunk number 8: Chapter05_LinearModel.prnw:210-211 ################################################### anova(model) ################################################### ### code chunk number 9: Chapter05_LinearModel.prnw:253-257 ################################################### summary(model)$r.squared R <- with(Orange,cor(circumference,age)) # R2(model,adj=FALSE) # not yet defined R^2 ################################################### ### code chunk number 10: Chapter05_LinearModel.prnw:319-320 (eval = FALSE) ################################################### ## levels(X) <- c("Afr","Ame","Asi") ################################################### ### code chunk number 11: Chapter05_LinearModel.prnw:329-330 (eval = FALSE) ################################################### ## X = factor(datasource, levels=c("never","often","always"), ordered=TRUE) ################################################### ### code chunk number 12: Chapter05_LinearModel.prnw:341-342 (eval = FALSE) ################################################### ## table(X) ################################################### ### code chunk number 13: Chapter05_LinearModel.prnw:358-359 (eval = FALSE) ################################################### ## contrasts(X)<-"contr.treatment" ################################################### ### code chunk number 14: Chapter05_LinearModel.prnw:401-405 ################################################### Strawberries = read.csv("Strawberries.csv") Strawberries X = acomp( Strawberries[,1:3] ) biomass = Strawberries$biomass ################################################### ### code chunk number 15: Chapter05_LinearModel.prnw:411-412 ################################################### Y = log(biomass) ################################################### ### code chunk number 16: Chapter05_LinearModel.prnw:444-452 ################################################### rescale = function(xin,rout,rin=range(xin, na.rm=FALSE)){ xout = rout[1] + (xin-rin[1])*(rout[2]-rout[1])/(rin[2]-rin[1]) return(xout) } Ycex = rescale(Y, c(0.5,2) ) Ycol = grey( rescale(Y,c(0.1,0.8)) ) plot(X,cex=Ycex) plot(X,pch=19,col=Ycol) ################################################### ### code chunk number 17: Chapter05_LinearModel.prnw:464-466 (eval = FALSE) ################################################### ## require("colorspace") ## Ycol = hex( mixcolor( rescale(Y,c(0,1)), RGB(1,0,0), RGB(0,0,1)) ) ################################################### ### code chunk number 18: XComp-VisualSizeColor ################################################### opar <- par(mar=c(3,3,0,0)) plot(X,pch=21,bg=Ycol,cex=Ycex) (Ylev <- pretty(Y)) legend(x=0,y=0.97,pch=21,legend=Ylev, pt.cex=rescale(Ylev, c(0.5,2),range(Y)), pt.bg=grey( rescale(Ylev,c(0.1,0.8),range(Y)) ) ) par(opar) ################################################### ### code chunk number 19: XComp-Visualpwlr ################################################### opar <- par(mfrow=c(3,3),mar=c(1,1,0.5,0.5), oma=c(3,3,0,0)) for(i in 1:3){ for(j in 1:3){ plot(log(X[,i]/X[,j]),Y,pch=ifelse(i!=j,19,"")) if(i==j){text(x=0,y=mean(range(Y)), labels=colnames(X)[i],cex=1.5)} } } mtext(text=c("pairwise log-ratio","dependent variable"),side=c(1,2),at=0.5,line=2,outer=TRUE) par(opar) ################################################### ### code chunk number 20: Chapter05_LinearModel.prnw:566-569 ################################################### (model = lm(Y~ilr(X))) (a = coef(model)[1]) (b = ilrInv(coef(model)[-1],orig=X)) ################################################### ### code chunk number 21: XComp-surfA ################################################### plot(X) straight(mean(X), b, lwd = 2, col = "black",lty=2) ################################################### ### code chunk number 22: XComp-surfB ################################################### myY = pretty(Y) refX = mean(X) + ((myY - a)/norm(b)^2) * b plot(refX, add = TRUE, pch = 19) ################################################### ### code chunk number 23: XComp-surfC ################################################### orthoComp = function(x) { ilrInv(ilr(x) %*% matrix(c(0, 1, -1, 0), ncol = 2)) } mygreyscale = grey((0:nrow(refX))/nrow(refX)) for (i in 1:nrow(refX)) { straight(acomp(refX[i, ]), orthoComp(b), col = mygreyscale[i],lwd = 2) } ################################################### ### code chunk number 24: XComp-surface ################################################### opar <- par(mar=c(3,3,0,1)) plot(X) straight(mean(X), b, lwd = 2, col = "black",lty=2) myY = pretty(Y) refX = mean(X) + ((myY - a)/norm(b)^2) * b plot(refX, add = TRUE, pch = 19) orthoComp = function(x) { ilrInv(ilr(x) %*% matrix(c(0, 1, -1, 0), ncol = 2)) } mygreyscale = grey((0:nrow(refX))/nrow(refX)) for (i in 1:nrow(refX)) { straight(acomp(refX[i, ]), orthoComp(b), col = mygreyscale[i],lwd = 2) } par(opar) ################################################### ### code chunk number 25: Chapter05_LinearModel.prnw:643-644 ################################################### (varb=ilrvar2clr(vcov(model)[-1,-1])) ################################################### ### code chunk number 26: ConfidenceX ################################################### par(mar=c(3,3,0,1)) names(b) = colnames(X) scaleB = 1 plot(scaleB*b) plot(0*b,add=TRUE,pch=20) alpha = 0.05 rF = sqrt(qf(1-alpha,nrow(varb)-1,model$df)) ellipses(scaleB*b,scaleB^2*varb,rF) ################################################### ### code chunk number 27: Chapter05_LinearModel.prnw:691-693 ################################################### newX <- acomp(c(1,2,3)) (newPredicted <- predict(model,newdata=list(X=newX))) ################################################### ### code chunk number 28: Chapter05_LinearModel.prnw:699-700 ################################################### predict(model, newdata = list(X = newX), interval = "confidence", level = 0.95) ################################################### ### code chunk number 29: Chapter05_LinearModel.prnw:705-710 ################################################### (newDesignmatrix = model.matrix(terms(model),data=list(X=newX,Y=0))) (errVar = newDesignmatrix %*% vcov(model)%*% t(newDesignmatrix)) alpha=0.05 (confidenceInterval = newPredicted + qt(c(alpha/2,1-alpha/2),model$df.residual)*sqrt(errVar)) (predictiveInterval = newPredicted + qt(c(alpha/2,1-alpha/2),model$df)*sqrt(summary(model)$sigma^2+errVar)) ################################################### ### code chunk number 30: plotPredVal ################################################### opar <- par(mar=c(3,3,0,0)) Predicted = predict(model) plot(Predicted,Y) abline(0,1) par(opar) ################################################### ### code chunk number 31: Chapter05_LinearModel.prnw:814-815 ################################################### cor(Y,Predicted,method="pearson") ################################################### ### code chunk number 32: Chapter05_LinearModel.prnw:819-821 ################################################### Residuals = model$residuals boxplot(list(Y-mean(Y),Predicted-mean(Y),Residuals)) ################################################### ### code chunk number 33: Chapter05_LinearModel.prnw:828-830 ################################################### summary(model) anova(model) ################################################### ### code chunk number 34: OutA ################################################### Predicted = predict(model) Residuals = resid(model) plot(Predicted,Residuals) ################################################### ### code chunk number 35: OutAll ################################################### par(mfrow=c(2,2)) plot(model,add.smooth=FALSE) ################################################### ### code chunk number 36: OutD ################################################### Leverages = hatvalues(model) boxplot(Leverages) CooksDistances = cooks.distance(model) plot(Predicted,CooksDistances) ################################################### ### code chunk number 37: Chapter05_LinearModel.prnw:1001-1006 ################################################### StudentizedResiduals = rstudent(model) qqnorm(StudentizedResiduals) qqline(StudentizedResiduals) shapiro.test(StudentizedResiduals) boxplot(StudentizedResiduals) ################################################### ### code chunk number 38: Chapter05_LinearModel.prnw:1027-1028 ################################################### cor.test(Predicted,abs(StudentizedResiduals),method="spearman") ################################################### ### code chunk number 39: Chapter05_LinearModel.prnw:1035-1038 ################################################### plot(Predicted, abs(StudentizedResiduals)) plot(Predicted, sqrt(abs(StudentizedResiduals))) plot(Predicted, log(abs(StudentizedResiduals))) ################################################### ### code chunk number 40: Chapter05_LinearModel.prnw:1049-1053 ################################################### require(robustbase) (modelRob = lmrob(Y~ilr(X))) (aRob = coef(modelRob)[1]) (bRob = ilrInv(coef(modelRob)[-1],orig=X)) ################################################### ### code chunk number 41: OutB ################################################### ResidualsRob= resid(modelRob) PredictedRob= Y-resid(modelRob) plot(PredictedRob,Residuals) ################################################### ### code chunk number 42: OutC ################################################### bxpstats <- boxplot(ResidualsRob,plot=FALSE) plot(PredictedRob,Y,pch=ifelse(ResidualsRob %in% bxpstats$out,33,20)) ################################################### ### code chunk number 43: Chapter05_LinearModel.prnw:1082-1083 ################################################### cor(Y,Predicted,method="spearman") ################################################### ### code chunk number 44: Chapter05_LinearModel.prnw:1151-1180 ################################################### sq <- function(x) { n <- if(length(dim(x))>1) ncol(x) else length(x) j <- rep(1:n,n) k <- rep(1:n,each=n) s <- j<=k f <- ifelse(j==k,1,2)[s] if(length(dim(x))>1) f*x[,j[s]]*x[,k[s]] else x[j[s]]*x[k[s]] } model2 = lm(Y~ilr(X)+sq(ilr(X))) model2 (a = coef(model2)[1]) (b = coef(model2)[2:ncol(X)]) unsq <- function(x) { #n(n-1)/2+n=l #n^2+n=2l #n^2+n-2l=0 #n=(-1+sqrt(1+8l))/2 n = (sqrt(8*length(x)+1)-1)/2 j <- rep(1:n,n) k <- rep(1:n,each=n) s <- j<=k i <- j+(k-1)*n C <- diag(n) C[i[s]]<-x C <- C+t(C) diag(C) <- diag(C)/2 C } (C = unsq(coef(model2)[-(1:ncol(X))])) ################################################### ### code chunk number 45: Chapter05_LinearModel.prnw:1196-1197 ################################################### (xOpt = ilrInv(solve(C,b)/2,orig=X)) ################################################### ### code chunk number 46: Chapter05_LinearModel.prnw:1201-1202 ################################################### eigen(C) ################################################### ### code chunk number 47: Chapter05_LinearModel.prnw:1211-1212 ################################################### anova(model2) ################################################### ### code chunk number 48: Chapter05_LinearModel.prnw:1238-1245 ################################################### #require(compositions) GeoChemSed=read.csv("GraVel.csv",header=TRUE) Y=acomp(GeoChemSed[,7:10]) names(GeoChemSed) Covariables=GeoChemSed[,c("discharge","relief","grain","position")] X1 = Covariables$discharge X2 = Covariables$relief ################################################### ### code chunk number 49: Chapter05_LinearModel.prnw:1264-1268 ################################################### X4 = factor(Covariables$position) X3 = factor(Covariables$grain,c("f","m","c"),ordered=TRUE) class(X3) X3 ################################################### ### code chunk number 50: Chapter05_LinearModel.prnw:1280-1283 ################################################### levels(X3) levels(X3) <- c("fine","medium","coarse") X3 ################################################### ### code chunk number 51: Chapter05_LinearModel.prnw:1289-1292 ################################################### contrasts(X3) contrasts(X3)<-"contr.treatment" contrasts(X3) ################################################### ### code chunk number 52: condVis1 ################################################### opar <- par(mar=c(4,4,1,1)) pairwisePlot(cbind(X1,X2,logX1=log(X1)),clr(Y)) par(opar) ################################################### ### code chunk number 53: condVis2 ################################################### opar = par(xpd=NA,no.readonly=TRUE) plot(Y,pch=c(3,20)[X4],col=c("red","green","blue")[X3]) legend(x=0.5,y=0.45,abbreviate(levels(X4), minlength=1), pch=c(3,20)) legend(x=0.75,y=0.0,abbreviate(levels(X3), minlength=1), pch=20,col=c("red","green","blue"),yjust=0) par(opar) ################################################### ### code chunk number 54: Chapter05_LinearModel.prnw:1392-1393 (eval = FALSE) ################################################### ## legend(locator(1),levels(X3),pch=20,col=c("red","green","blue"),xpd=NA,yjust=0) ################################################### ### code chunk number 55: Chapter05_LinearModel.prnw:1401-1402 ################################################### plot(Y,pch=as.numeric(X4),col=X3) ################################################### ### code chunk number 56: boxmulti ################################################### boxplot(Y,X3,notch=TRUE) ################################################### ### code chunk number 57: fullPairwise ################################################### opar <- par(mar=c(4,4,0,0)) pairwisePlot(cbind(X1,X2,logX2=log(X2)),clr(Y),pch=c(3,20)[X4],col=c("red","green","blue")[X3]) par(opar) ################################################### ### code chunk number 58: Chapter05_LinearModel.prnw:1491-1492 ################################################### (model = lm(ilr(Y)~log(X1) )) ################################################### ### code chunk number 59: Chapter05_LinearModel.prnw:1498-1499 ################################################### ilrInv(coef(model),orig=Y) ################################################### ### code chunk number 60: Chapter05_LinearModel.prnw:1509-1510 ################################################### anova(model) ################################################### ### code chunk number 61: Chapter05_LinearModel.prnw:1516-1518 ################################################### model=lm(ilr(Y)~log(X2)) anova(model) ################################################### ### code chunk number 62: Chapter05_LinearModel.prnw:1562-1568 ################################################### contrasts(X3) <- "contr.treatment" (model = lm(ilr(Y)~X3)) (a = ilrInv(coef(model)[1,],orig=Y)) (b = ilrInv(rbind(0,coef(model)[-1,]),orig=Y)) (mu = a+b) (Sigma = ilrvar2clr(var(model)) ) ################################################### ### code chunk number 63: Chapter05_LinearModel.prnw:1574-1577 ################################################### plot(mu,pch=20) plot(Y,pch=".",add=TRUE) ellipses(mu,Sigma,2) ################################################### ### code chunk number 64: Chapter05_LinearModel.prnw:1587-1588 ################################################### anova(model) ################################################### ### code chunk number 65: Chapter05_LinearModel.prnw:1606-1607 ################################################### (model = lm(ilr(Y)~log(X1)+X2+X3+X4 )) ################################################### ### code chunk number 66: Chapter05_LinearModel.prnw:1632-1633 ################################################### coef(model) ################################################### ### code chunk number 67: Chapter05_LinearModel.prnw:1637-1638 ################################################### (coefs <- ilrInv(coef(model),orig=Y)) ################################################### ### code chunk number 68: dispCoef ################################################### barplot(coefs,las=2,xlim=c(0,11)) ################################################### ### code chunk number 69: Chapter05_LinearModel.prnw:1654-1655 (eval = FALSE) ################################################### ## vcov(model) ################################################### ### code chunk number 70: Chapter05_LinearModel.prnw:1663-1665 ################################################### cat(paste(sapply(1:nrow(coefs),function(i){paste(i, rownames(coefs)[i],sep=" = ")}), collapse="; ")) ################################################### ### code chunk number 71: coefell ################################################### #vcovAcomp <- function(object,...) { # co <- coef(object) # aperm(structure(vcov(object,...),dim=c(dim(co),dim(co)),dimnames=c(dimnames(co),dimnames(co))),c(2,4,1,3)) #} #qHotellingsT2 <- function(x,p,m) qf(x,p,m-p+1)*(p*m)/(m-p+1) #ConfRadius <- function(model,prob=1-alpha,alpha) { # sqrt(qHotellingsT2(prob,ncol(coef(model)),model$df.residual)) #} vars <- vcovAcomp(model) dim(vars) alpha=0.05 plot(coefs,pch=as.character(1:nrow(coefs))) plot(coefs-coefs,pch=20,add=TRUE) for(i in 1:nrow(coefs)){ ellipses(acomp(coefs[i,,drop=FALSE]), ilrvar2clr(vars[,,i,i]), r=ConfRadius(model,1-alpha), col="gray") } par(xpd=FALSE) ################################################### ### code chunk number 72: lm-biplot ################################################### B = clr(coefs[-1,]) svdlm = svd(B) opar <- par(xpd=NA,mar=c(2,2,0,0)) coloredBiplot(svdlm, scale =0, xarrows=TRUE, ypch=4, ynames=colnames(B), xnames=rownames(B)) sum(svdlm$d[1:2]^2)/sum(svdlm$d^2) par(opar) ################################################### ### code chunk number 73: Chapter05_LinearModel.prnw:1805-1806 ################################################### anova(model) ################################################### ### code chunk number 74: Chapter05_LinearModel.prnw:1889-1893 ################################################### anova(lm(ilr(Y)~X2 + X3 + X4 + log(X1))) [5,] anova(lm(ilr(Y) ~ log(X1) + X3 + X4 + X2)) [5,] anova(lm(ilr(Y) ~ log(X1) + X2 + X4 + X3)) [5,] anova(lm(ilr(Y) ~ log(X1) + X2 + X3 + X4)) [5,] ################################################### ### code chunk number 75: Chapter05_LinearModel.prnw:1915-1916 ################################################### Pred = ilrInv(predict(model),orig=Y) ################################################### ### code chunk number 76: CompPred (eval = FALSE) ################################################### ## plot(Y,pch=".") ## plot(Predicted,add=TRUE,pch=20) ## segments(x0=Predicted,y=Y) ################################################### ### code chunk number 77: predobs ################################################### #opar <- par(oma=c(0,0,2,0),mar=c(4,4,1,0)) #pairwisePlot(clr(ilrInv(predict(model),orig=Y)),clr(Y)) #title("Observed vs. Predicted",outer=TRUE) #par(opar) opar <- par(mfrow=c(4,4), mar=c(2,2,1,1), oma=c(4,4,0,0)) Ypred = ilrInv(predict(model),orig=Y) for(i in 1:4){ for(j in 1:4){ plot(log(Y[,i]/Y[,j]), log(Ypred[,i]/Ypred[,j]), pch=ifelse(i!=j,1,"")) if(i==j){text(x=0,y=0, labels=colnames(Y)[i],cex=1.5) }else{ abline(a=0,b=1,col="gray",lwd=3) } } } mtext(text=c("observed log-ratio","predicted log-ratio"),side=c(1,2),at=0.5,line=2,outer=TRUE) par(opar) ################################################### ### code chunk number 78: Chapter05_LinearModel.prnw:1996-2000 ################################################### newdata=list(X1=geometricmean(X1), X2=mean(X2), X3=factor("coarse"), X4=factor("north")) ################################################### ### code chunk number 79: Chapter05_LinearModel.prnw:2008-2009 ################################################### prediction <- ilrInv(predict(model,newdata=newdata),orig=Y) ################################################### ### code chunk number 80: Chapter05_LinearModel.prnw:2088-2089 ################################################### (varEpsilon = var(model)) ################################################### ### code chunk number 81: Chapter05_LinearModel.prnw:2133-2146 ################################################### getModelMatrix <- function(object,newdata=NULL,na.action=na.pass) { if( is.null(newdata) ) return(model.matrix(object)) Terms <- delete.response(terms(object)) mf <- model.frame(Terms,newdata,na.action=na.action,xlev = object$xlevels) if (!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, mf) model.matrix(Terms, mf, contrasts.arg = object$contrasts) } (X <- getModelMatrix(model,newdata)) XX <- kronecker(diag(ncol(predict(model))),X) (estVar = XX %*% vcov(model) %*% t(XX)) (predVar = estVar+varEpsilon) ################################################### ### code chunk number 82: apredCI ################################################### alpha=0.05 plot(Y,pch=".") plot(prediction,pch=20,add=TRUE) ellipses(prediction,ilrvar2clr(estVar), r=ConfRadius(model,1-alpha)) ellipses(prediction,ilrvar2clr(predVar), r=ConfRadius(model,1-alpha)) #ellipses(prediction,ilrvar2clr(varEpsilon),r=ConfRadius(model,1-alpha),col="blue") #title("Prediction CI and Predictive Interval",outer=TRUE) ################################################### ### code chunk number 83: Chapter05_LinearModel.prnw:2197-2198 ################################################### Resid = ilrInv(resid(model),orig=Y) ################################################### ### code chunk number 84: residTer ################################################### plot(Resid, cex=0.5) #title("Ternary Diagrams of Residuals",outer=TRUE,line=-1) ################################################### ### code chunk number 85: residTer2 (eval = FALSE) ################################################### ## plot(4*Resid) ################################################### ### code chunk number 86: resbox ################################################### boxplot(Resid) #title("Boxplot of Residuals",outer=TRUE,line=-1) ################################################### ### code chunk number 87: residQQnorm ################################################### qqnorm(Resid) #title("Normal QQ-Plot of Residuals",outer=TRUE,line=-1) ################################################### ### code chunk number 88: respred ################################################### opar <- par(oma=c(3,3,0,0),mar=c(4,4,1,0)) pairwisePlot(clr(Pred),clr(Resid)) mtext(text=c("predicted values (clr)","residuals (clr)"),side=c(1,2),at=0.5,line=2,outer=TRUE) #title("Residuals vs. Predicted",outer=TRUE) par(opar) ################################################### ### code chunk number 89: respred2 ################################################### opar <- par(mfrow=c(4,4),mar=c(2,2,1,1), oma=c(4,4,0,0)) for(i in 1:4){ for(j in 1:4){ plot(log(Pred[,i]/Pred[,j]),log(Resid[,i]/Resid[,j]) ,pch=ifelse(i!=j,19,"")) if(i==j){text(x=0,y=0,labels=colnames(Resid)[i],cex=1.5)} else{abline(h=0)} } } mtext(text=c("predicted values","residuals"),side=c(1,2),at=0.5,line=2,outer=TRUE) par(opar) ################################################### ### code chunk number 90: Chapter05_LinearModel.prnw:2397-2401 (eval = FALSE) ################################################### ## opar <- par(oma=c(0,0,2,0),mar=c(4,4,1,0)) ## pairwisePlot(clr(Pred),clr(Resid),col=as.numeric(X3)) ## legend(locator(1),levels(X3),col=as.numeric(1:3),pch=1,xpd=NA) ## par(opar) ################################################### ### code chunk number 91: Chapter05_LinearModel.prnw:2456-2470 ################################################### #R2 <- function(object,...,adj=TRUE,ref=0) { # if( !is.numeric(ref)) # ref<-Recall(ref,...,adj=adj) # pr <- predict(object) # re <- resid(object) # dfres<-object$df.residual # n <- nrow(re) # y <- pr+re # erg <- if( adj ) # 1-(mvar(re)*(n-1)/dfres)/mvar(y) # else # 1-mvar(re)/mvar(y) # (erg-ref)/(1-ref) #} ################################################### ### code chunk number 92: Chapter05_LinearModel.prnw:2480-2489 ################################################### #R2(model) R2(lm(ilr(Y)~log(X1)+X2)) R2(lm(ilr(Y)~log(X1)+X3)) R2(lm(ilr(Y)~log(X1)+X2+X4)) R2(lm(ilr(Y)~log(X1)+X2+X3)) R2(lm(ilr(Y)~X2+X3+X4)) R2(lm(ilr(Y)~log(X1)+X2+X3+X4)) #R2(lm(ilr(Y)~log(X1)*X3)) #R2(lm(ilr(Y)~log(X1)*X2*X3*X4)) ################################################### ### code chunk number 93: Chapter05_LinearModel.prnw:2532-2536 ################################################### compY <- acomp(Y[GeoChemSed$grain=="c",]) compX <- acomp(Y[GeoChemSed$grain=="m",]) dim(compY) dim(compX) ################################################### ### code chunk number 94: compcomp ################################################### opar<-par(mar=c(4,4,0,0),oma=c(3,3,0.1,0.1)) pairwisePlot(clr(compX),clr(compY)) mtext(text=c("medium-sized","coarse"),side=c(1,2),at=0.5,line=2,outer=TRUE) par(opar) #par(mfrow=c(4,4),mar=c(2,2,1,1), oma=c(4,4,0,0)) #for(i in 1:4){ # for(j in 1:4){ # plot(log(compX[,i]/compX[,j]),log(compY[,i]/compY[,j]) ,pch=ifelse(i!=j,19,"")) # if(i==j){text(x=0,y=0,labels=colnames(compX)[i],cex=1.5)} # } #} #mtext(text=c("medium-sized","coarse"),side=c(1,2),at=0.5,line=2,outer=TRUE) ################################################### ### code chunk number 95: Chapter05_LinearModel.prnw:2580-2582 ################################################### (modelCC <- lm(ilr(compY)~ilr(compX))) anova(modelCC) ################################################### ### code chunk number 96: Chapter05_LinearModel.prnw:2589-2590 ################################################### ilrInv(coef(modelCC)[1,],orig=compY) ################################################### ### code chunk number 97: biplot-selflm ################################################### (Bsvd <- svd( ilrvar2clr(coef(modelCC)[-1,]) )) opar <- par(xpd=NA,mar=c(2,2,0,1)) coloredBiplot(Bsvd,xarrows=TRUE,xnames=colnames(Y),ynames=colnames(Y)) par(opar) ################################################### ### code chunk number 98: Chapter05_LinearModel.prnw:2620-2625 ################################################### u = clrInv(t(Bsvd$u)) v = clrInv(t(Bsvd$v)) colnames(u) <- colnames(v) <- colnames(Y) print(u) print(v) ################################################### ### code chunk number 99: Chapter05_LinearModel.prnw:2694-2695 ################################################### summary(modelCC) ################################################### ### code chunk number 100: Chapter05_LinearModel.prnw:2717-2719 ################################################### (coorindex = factor(rep(1:ncol(ilr(compX)),each=nrow(compX)))) (simplemodel = lm( c(ilr(compY)) ~ c(ilr(compX)) + coorindex ) ) ################################################### ### code chunk number 101: Chapter05_LinearModel.prnw:2838-2850 ################################################### GeoChemSed=read.csv("GraVel.csv",header=TRUE) Y=acomp(GeoChemSed[,7:10]) names(GeoChemSed) Covariables=GeoChemSed[,c("discharge","relief","grain","position")] X1 = Covariables$discharge X2 = Covariables$relief X4 = factor(Covariables$position) X3 = factor(Covariables$grain,c("f","m","c"),ordered=TRUE) levels(X3) <- c("fine","medium","coarse") contrasts(X3)<-"contr.treatment" anova(lm(ilr(Y)~log(X1)+ X2 + X3 + X4 + log(X1)*X2+X2*X3+X3*X4)) ################################################### ### code chunk number 102: Chapter05_LinearModel.prnw:2857-2858 (eval = FALSE) ################################################### ## anova(lm(ilr(Y)~log(X1)*X2+X2*X3+X3*X4)) ################################################### ### code chunk number 103: Chapter05_LinearModel.prnw:2862-2863 (eval = FALSE) ################################################### ## anova(lm(ilr(Y)~log(X1)*X2+X3*X4)) ################################################### ### code chunk number 104: Chapter05_LinearModel.prnw:2866-2867 ################################################### anova(lm(ilr(Y)~log(X1)*X2+X2*X4+log(X1)*X3*X4)) ################################################### ### code chunk number 105: Chapter05_LinearModel.prnw:2875-2876 ################################################### anova(lm(ilr(Y)~log(X1)*X2+X2*X4+log(X1):X3:X4)) ################################################### ### code chunk number 106: Chapter05_LinearModel.prnw:2946-2952 ################################################### AIC(model) AIC(lm(ilr(Y)~log(X1)+X3)) AIC(lm(ilr(Y)~log(X1)+X2+X3+X4)) AIC(lm(ilr(Y)~log(X1)*X3)) AIC(lm(ilr(Y)~log(X1)*X2*X3*X4)) AIC(lm(ilr(Y)~log(X1)*X2*X3+X1*X3*X4)) ################################################### ### code chunk number 107: Chapter05_LinearModel.prnw:3042-3044 ################################################### maxmodel = lm(ilr(Y)~log(X1)*X3+log(X1)*X4+X3*X4 ) (terms = rownames(coef(maxmodel))[-1]) ################################################### ### code chunk number 108: Chapter05_LinearModel.prnw:3048-3051 ################################################### X3medium = X3=="medium" X3coarse = X3=="coarse" X4south = X4=="south" ################################################### ### code chunk number 109: Chapter05_LinearModel.prnw:3056-3057 ################################################### idx = sapply(1:length(terms), function(i){combn(terms,i) } ) ################################################### ### code chunk number 110: Chapter05_LinearModel.prnw:3070-3080 ################################################### res = lapply(idx, function(termstaken){ tt = termstaken if(is.null(dim(tt)))dim(tt)=c(length(tt),1) lapply(1:ncol(tt), function(icol){ frm = paste( tt[,icol], collapse="+") frm = paste("ilr(Y)~", frm, sep="") mylm = lm(as.formula(frm)) return(mylm) }) }) ################################################### ### code chunk number 111: Chapter05_LinearModel.prnw:3086-3091 ################################################### length(res) length(res[[1]]) res = unlist(res, recursive=FALSE) length(res) names(res) = 1:length(res) ################################################### ### code chunk number 112: Chapter05_LinearModel.prnw:3103-3121 ################################################### mylogLik.mlm = function(mylm){ x = mylm$residuals detvarx = ifelse(length(dim(x))>1, determinant(var(x))$modulus, log(var(x))) ncolx = ifelse(length(dim(x))>1,ncol(x),1) -0.5*(sum(mahalanobis(x,0*mean(x),var(x)))+detvarx + ncolx * log(2*pi)) } #( lL = sapply(res,logLik) ) lL = sapply(res,mylogLik.mlm) lc = sapply(res,function(mylm){length( coef(mylm)) }) rvc = sapply(res,function(mylm){ d = max(ncol(mylm$residuals),1) #$ return(d*(d+1)/2) }) np = lc + rvc sort(myAIC <- -2*lL + 2*np)[1:10] sort(myBIC <- -2*lL + log(nrow(Y)) * np)[1:10] ################################################### ### code chunk number 113: Chapter05_LinearModel.prnw:3125-3126 ################################################### lapply(res[c(3,4,2,1)], coef) ################################################### ### code chunk number 114: Chapter05_LinearModel.prnw:3139-3140 ################################################### maxmodel = lm( ilr(Y)~log(X1)+X2+X3+X4) ################################################### ### code chunk number 115: Chapter05_LinearModel.prnw:3144-3145 ################################################### (idx = lapply(2:(ncol(Y)-1), function(i){combn(ncol(Y),i)})) ################################################### ### code chunk number 116: Chapter05_LinearModel.prnw:3153-3165 ################################################### res = lapply(idx, function(subcomps){ lapply(1:ncol(subcomps), function(icol){ idparts = unique(c(subcomps[,icol],1:ncol(Y))) sY = acomp(Y,idparts) Y = ilr(sY) Y1 = Y[,1:(nrow(subcomps)-1)] Y2 = Y[,-(1:(nrow(subcomps)-1))] lm1 = lm( Y1~log(X1)+X2+X3+X4) lm2 = lm( Y2~1 ) return( list(lm=lm1,nolm=lm2) ) }) }) ################################################### ### code chunk number 117: Chapter05_LinearModel.prnw:3172-3182 ################################################### length(res) length(res[[1]]) res = unlist(res, recursive=FALSE) length(res) names(res) = unlist( lapply(idx,function(subcomps){ lapply(1:ncol(subcomps), function(icol){ paste(colnames(Y)[subcomps[,icol]],collapse="-") } ) } ) ) ################################################### ### code chunk number 118: Chapter05_LinearModel.prnw:3191-3209 ################################################### mylogLik = function(mylm){ x = cbind(mylm$lm$residuals,mylm$nolm$residuals) detvarx = ifelse(length(dim(x))>1, determinant(var(x))$modulus, log(var(x))) ncolx = ifelse(length(dim(x))>1,ncol(x),1) -0.5*(sum(mahalanobis(x,0*mean(x),var(x)))+detvarx + ncolx * log(2*pi)) } myPars = function(mylm){ x = cbind(mylm$lm$residuals,mylm$nolm$residuals) d = ncol(x) return(d*(d+1)/2 + length( coef(mylm$lm)) + length( coef(mylm$nolm) )) } ( lL = sapply(res,mylogLik) ) ( np = sapply(res,myPars) ) sort(myAIC <- -2*lL + 2*np) sort(myBIC <- -2*lL + log(nrow(Y)) * np) ################################################### ### code chunk number 119: Chapter05_LinearModel.prnw:3217-3220 ################################################### incrlL = mylogLik.mlm(maxmodel) - lL incrnp = ( length(coef(maxmodel)) + 3*4/2 ) - np (alpha.levels <- pchisq(-2*incrlL,df=incrnp*2,lower.tail=FALSE))