An Exploration of Virtual Auditory Shape Perception

Appendix A:
Regression Algorithm

	I performed the analysis of the screening test data using the S Plus software.
In order to handle the vertical confusions in the height estimation data, I
need to modify some of the S Plus analysis routines, and write some of my own
functions.  This code appears below.


#A simple function to compute the residues given the slope

# and intercept of a line.

residuate <- function(stimuli, response, b, m) {(response - (m * stimuli +
b))}


# This is a weighting function for robust regression that 

# discounts points that are closer to a line with the negative

# slope of the previous iteration of the regression.

wt.signerr <-function(x,y,coef)

{

	fr <-  (residuate(x,y, coef[1], coef[2]))

	tvr <- (residuate(x,y,coef[1], (-1*coef[2])))

	w<- ifelse( (abs(fr) < abs(tvr)), 1.0, 0.0)

	w

}


# This is a modified version of SPlus' Robust Regression Function.

# It employs the weighting function for vertical confusions.

rxreg <- function(x, y, w = rep(1, nrow(x)), int = TRUE, init = lsfit(x, y,
w, int = 

	FALSE)$coef, method = wt.signerr, wx, iter = 20, acc = 10 * .Machine$

	single.eps^0.5, test.vec = "resid")

{

	   objctv<-x

	   subjctv<-y

	irls.delta <- function(old, new)

	{

		a <- sum((old - new)^2)

		b <- sum(old^2)

		if(b >= 1 || a < b * .Machine$double.xmax)

			sqrt(a/b)

		else .Machine$double.xmax

	}

	irls.rrxwr <- function(x, w, r)

	{

		w <- sqrt(w)

		max(abs((as.vector(r * w) %*% x)/sqrt(as.vector(w) %*% (x^2))))/

			sqrt(sum(w * r^2))

	}

	if(!(any(test.vec == c("resid", "coef", "w", "NULL")) || is.null(

		test.vec)))

		stop("invalid testvec")

	if(int)

		x <- cbind(1, x)

	else x <- as.matrix(x)

	if(!missing(wx)) {

		if(length(wx) != nrow(x))

			stop("Length of wx must equal number of observations")

		if(any(wx < 0))

			stop("Negative wx value")

		w <- w * wx

	}

	coef <- init

	if(ncol(x) != length(coef))

		stop("Must have same number of initial values as coefficients")

	resid <- y - x %*% coef

	converged <- FALSE

	status <- "converged"

	conv <- NULL

	method.in.control <- method.exit <- FALSE

	for(iiter in 1:iter) {

		if(!is.null(test.vec))

			previous <- get(test.vec)

		scale <- median(abs(resid))/0.6745

		if(scale == 0) {

			convi <- 0

			method.exit <- TRUE

			status <- "could not compute scale of residuals"

		}

		else {

			w <- method(objctv,subjctv,coef)

			if(!missing(wx))

				w <- w * wx

			temp <- lsfit(x, y, w, int = FALSE)

			coef <- temp$coef

			resid <- temp$residuals

			if(!is.null(test.vec))

				convi <- irls.delta(previous, get(test.vec))

			else convi <- irls.rrxwr(x, w, resid)

		}

		conv <- c(conv, convi)

		converged <- convi <= acc

		done <- method.exit || (converged && !method.in.control)

		if(done)

			break

	}

	if(!done)

		warning(status <- paste("failed to converge in", iter, "steps")

			)

	if(!missing(wx)) {

		tmp <- (wx != 0)

		w[tmp] <- w[tmp]/wx[tmp]

	}

		 list(coef = coef, residuals = resid, w = w, int = int, conv = conv, 

		status = status)

}


# Returns robust regression of data.

confusion <- function(name)

{

	regval <-  rxreg(height(name), estimate(name))

	text(height(name)[regval$w==1], estimate(name)[regval$w==1],"*")

	text(height(name)[regval$w==0], estimate(name)[regval$w==0], "0")

	abline(regval)

	abline(regval$coef[1],(-1*regval$coef[2]))

	axes("Confusion Plot of Height Estimation", paste("Subject:", name), "Height",
"Estimate")

	regval

}


# Plots the robust regression of data outside a core area

# of expected heights.  Also returns regression information.

cconfusion <- function(name,core)

{

	ht<-height(name)[abs(height(name))>core]

	est<-estimate(name)[abs(height(name))>core]

	regval <-  rxreg(ht, est)


	axes("Confusion Plot of Height Estimation", paste("Subject:", name), "Height",
"Estimate", xlim=c(-13,13), ylim=c(-13,13))

	text(ht[regval$w==1], est[regval$w==1],"X")

	text(ht[regval$w==0], est[regval$w==0], "0")

	abline(regval)

	abline(regval$coef[1],(-1*regval$coef[2]))

	    regval

}


#Plots residues of regressed data, resolving vertical confusions.

cconfresplot <- function(name,core)

{

	ht<-height(name)[abs(height(name))>core]

	est<-estimate(name)[abs(height(name))>core]

	regval <-  rxreg(ht, est)


	est[regval$w==0]<- -1*(est[regval$w==0]-regval$coef[1])

	leftovers<-residuate(ht, est, regval$coef[1], regval$coef[2])

	plot(ht, leftovers, main="Residue of Height Estimate with Confusions
Resolved", sub=paste("Subject:", name), xlab="Height", ylab="Residue")


}


#Returns residues of regressed data, resolving vertical confusions.

cconfres <- function(name,core)

{

	ht<-height(name)[abs(height(name))>core]

	est<-estimate(name)[abs(height(name))>core]

	regval <-  rxreg(ht, est)


	est[regval$w==0]<- -1*(est[regval$w==0]-regval$coef[1])

	leftovers<-residuate(ht, est, regval$coef[1], regval$coef[2])

	leftovers

}


#Generates a log-log plot of Height Estimates

# with vertical confusions noted.


cconflog <- function(name,core)

{

	ht<-height(name)[abs(height(name))>core]

	est<-estimate(name)[abs(height(name))>core]

	regval <-  rxreg(ht, est)


	axes("log-log of Confusion Plot of Height Estimation", paste("Subject:",
name), "Height", "Estimate", xlim=c(-2,2), ylim=c(-2,2))

	text(logify(ht[regval$w==1]),logify( est[regval$w==1]),"X")

	text(logify(ht[regval$w==0]), logify(-1*(est[regval$w==0]-regval$coef[1])),
"0")

	est[regval$w==0]<- -1*(est[regval$w==0]-regval$coef[1])

	est

}


# This function returns a cored confusion count:

# A count of vertical confusions for data points.

ccc <- function(name,core) 

{

	ht<-height(name)[abs(height(name))>core]

	est<-estimate(name)[abs(height(name))>core]

	regval <-  rxreg(ht, est)

	length(regval$w[regval$w==0])

}