--- title: "MLR with delta mehod" author: "M. Shotwell" date: "February 2, 2018" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ### Function to compute finite difference gradients ```{r} ## compute gradient of probability function for a sequence of c1 at median c2 ## compute finite-difference gradient (c.f., nlme::fdHess) fdGrad <- function (pars, fun, ..., .relStep = (.Machine$double.eps)^(1/2), minAbsPar = 0) { ##pars <- as.numeric(pars) npar <- length(pars) incr <- ifelse(abs(pars) <= minAbsPar, .relStep, (abs(pars)-minAbsPar) * .relStep) ival <- do.call(fun, list(pars, ...)) diff <- rep(0,npar) sapply(1:npar, function(i) { del <- rep(0,npar) del[i] <- incr[i] (do.call(fun, list(pars+del, ...))-ival)/incr[i] }) } ``` ### Compue canonical variables from RRLDA ```{r} library('ElemStatLearn') library('nnet') ## read vowel training and testing data from HTF website data(vowel.train); data(vowel.test) dat <- vowel.train tst <- vowel.test ## compute canonical variables from reduced rank LDA ## do logistic regression with first two canonidal variables ## names of input variables inp <- paste("x", 1:10, sep=".") ## get input matrix inpx <- as.matrix(dat[inp]) ## center inputs globally inpx <- t(t(inpx) - colMeans(inpx)) ## compute within-cluster means cenx <- t(sapply(unique(dat$y), function(cls) colMeans(inpx[dat$y==cls,]))) ## compute pooled within-cluster var-cov sigx <- Reduce(`+`, lapply(unique(dat$y), function(cls) { x <- inpx[dat$y==cls,] x <- t(t(x)-colMeans(x)) t(x)%*%x }))/(nrow(dat)-11) ## 'sphere' the inputs and centers svdx <- svd(sigx) sphx <- inpx%*%svdx$v%*%diag(sqrt(1/svdx$d)) sphc <- cenx%*%svdx$v%*%diag(sqrt(1/svdx$d)) ## compute canonical inputs and centers svdc <- svd(sphc) canx <- sphx%*%svdc$v canc <- sphc%*%svdc$v canx_df <- data.frame(canx) colnames(canx_df) <- paste0('c',1:10) dat <- cbind(dat, canx_df) ``` ### Plot probability of class 1 membership using delta method interval ```{r} ## fit (multinomial) logistic regression fit <- multinom(y ~ c1 + c2, data=dat, Hess=TRUE) summary(fit) ## compute model matrix ## compute linear predictors for class membership linpred <- function(beta, x) { xb <- x %*% t(beta) return(xb) } ## compute class membership probabilities from linear predictors sigmoid <- function(xb) { pr <- exp(xb)/(1+rowSums(exp(xb))) pr <- cbind('1'=1-rowSums(pr), pr) return(pr) } ## compute approximate variance-covariance of \hat{\beta} vcv_beta <- solve(fit$Hessian) ## generate sequence of values of c1 at 5th percentile of c2 pre_dat <- data.frame(c1=seq(min(dat$c1), max(dat$c1), length.out=200), c2=quantile(dat$c2, 0.05)) ## compute associated design matrix pre_mmx <- model.matrix(~c1+c2, pre_dat) ## show region of interest plot(dat$c1, dat$c2, col=dat$y, xlab="Canonical Variable 1", ylab="Canonical Variable 2") abline(h=quantile(dat$c2, 0.05), lwd=2, lty=2) ## plot probabilities plot(pre_dat$c1, sigmoid(linpred(coef(fit),pre_mmx))[,1], xlab='Canonical Variable 1', ylab='Probability of Class 1', type='l', lwd=2) ## compute FD gradient of probability function at estimate of beta grd_pr1 <- fdGrad(coef(fit), function(beta) { sigmoid(linpred(beta, pre_mmx))[,1] }) ## compute approximate variance covariance of probabilities vcv_pr1 <- grd_pr1%*%vcv_beta%*%t(grd_pr1) ## compute probabilities with approximate 95% CI using delta method pre_dat$pr1_est <- sigmoid(linpred(coef(fit),pre_mmx))[,1] pre_dat$pr1_lo95 <- pmax(0, pre_dat$pr1_est - qnorm(0.975)*sqrt(diag(vcv_pr1))) pre_dat$pr1_hi95 <- pmin(1, pre_dat$pr1_est + qnorm(0.975)*sqrt(diag(vcv_pr1))) ## plot probabilities plot(pre_dat$c1, pre_dat$pr1_est, ylim=range(c(pre_dat$pr1_lo95,pre_dat$pr1_hi95)), xlab='Canonical Variable 1', ylab='Probability of Class 1', type='l', lwd=2) polygon(c(pre_dat$c1, rev(pre_dat$c1)), c(pre_dat$pr1_hi95, rev(pre_dat$pr1_lo95)),density = 15) ``` ### Plot probability of class 1 membership with bootstraps ```{r} ## function compute probability of class 1 membership prob_1 <- function(fit, dat) predict(fit, newdata=dat, type='prob')[,1] ## function to resample rows of data resample <- function(dat) { idx <- sample(1:nrow(dat), size=nrow(dat), replace=T) return(dat[idx,]) } pre_boot <- replicate(100, { dat <- resample(dat) fit <- multinom(y ~ c1 + c2, data=dat, trace=FALSE) prob_1(fit, pre_dat) }) ## plot probabilities plot(pre_dat$c1, prob_1(fit, pre_dat), ylim=range(pre_boot), xlab='Canonical Variable 1', ylab='Probability of Class 1', type='l', lwd=2) for(i in 1:ncol(pre_boot)) lines(pre_dat$c1, pre_boot[,i], col='#00000033') ```