Variable Signal
Nina Zumel
July 31, 2015
Functions
library(ggplot2)
# get the theoretical significance glm model (wrt deviance)
get_significance = function(model) {
delta_deviance = model$null.deviance - model$deviance
df = model$df.null - model$df.residual
sig = pchisq(delta_deviance, df, lower.tail=FALSE)
}
# return a frame of the deviance scores on the permuted data
permutation_test = function(dataf, ycol, nperm) {
nrows = dim(dataf)[1]
y = dataf[[ycol]]
X = dataf[, setdiff(colnames(dataf), ycol), drop=FALSE]
varnames = colnames(X)
fmla = paste("y", paste(varnames, collapse=" + "), sep=" ~ ")
deviances = numeric(nperm)
for(i in seq_len(nperm)) {
# random order of rows
ord = sample.int(nrows, size=nrows, replace=FALSE)
model = glm(fmla, data=cbind(y=y[ord], X),
family=binomial(link="logit"))
#print(summary(model))
deviances[[i]] =model$deviance
}
deviances
}
score_variable = function(dframe, ycol, var, nperm,
title='') {
df=data.frame(y=dframe[[ycol]], x=dframe[[var]])
mod = glm("y~x", data=df,
family=binomial(link="logit"))
vdev = mod$deviance
vperm = permutation_test(df, "y", nperm)
print(summary(mod))
# count how many times vdev >= deviances from perm test
num = sum(vperm <= vdev)
vscore = num/nperm
print(ggplot(data.frame(nullperm=vperm), aes(x=nullperm)) +
geom_density() + geom_vline(xintercept=vdev, color='red') +
ggtitle(paste(title, "left tail area ~", vscore)))
print(paste("Chi-sq estimate:", get_significance(mod)))
}
# ----
# for making data sets with both noise and weakly correlated variables
# ----
mkCoefs = function(ngood) {
nGoodN = ceiling(ngood/2)
nGoodC = ngood-nGoodN
coefs = list()
if(nGoodN>0) {
cN = lapply(seq_len(nGoodN),function(i) {
v = list()
v[[paste('gn', i, sep='_')]] = rnorm(1)
v
})
coefs = unlist(cN,recursive=FALSE)
}
if(nGoodN>0) {
cC = lapply(seq_len(nGoodC),function(i) {
v = list()
v[[paste('gc', i, sep='_')]] = list('a'=rnorm(1),
'b'=rnorm(1),
'c'=rnorm(1))
v
})
coefs = append(coefs,unlist(cC,recursive=FALSE))
}
coefs
}
# generate columns
genColumns = function(nrow,prefix,nNum,nCat) {
cols = list()
if(nNum>0) {
numCols = lapply(seq_len(nNum),function(i) rnorm(nrow))
names(numCols) = paste(prefix,'n_',seq_len(nNum),sep='')
cols = append(cols,numCols)
}
if(nCat>0) {
cCols = lapply(seq_len(nCat),function(i) {
sample(c('a','b','c'),nrow,replace=TRUE)
})
names(cCols) = paste(prefix,'c_',seq_len(nCat),sep='')
cols = append(cols,cCols)
}
data.frame(cols,stringsAsFactors=FALSE)
}
# evaluate coefs on a data frame
applyCoefs = function(coefs,d) {
v = numeric(nrow(d))
for(n in names(coefs)) {
cf = coefs[[n]]
if(is.list(cf)) {
# categorical
v = v + as.numeric(cf[d[,n,drop=TRUE]])
} else {
# numeric
v = v + cf[[1]]*d[,n,drop=TRUE]
}
}
v
}
# build a data frame with pure noise columns
# and columns weakly correlated with y
mkData = function(nrows, coefs, nnoise) {
noiseMagnitude = 1
d = data.frame(y = noiseMagnitude*rnorm(nrows),
stringsAsFactors = FALSE)
ngood = length(coefs)
if(ngood>0) {
ngC = sum(vapply(coefs,is.list,numeric(1)))
ngN = ngood - ngC
gd = genColumns(nrows,'g',ngN,ngC)
d = cbind(d,gd)
d$y = d$y + applyCoefs(coefs,d)
}
if(nnoise > 0) {
nnN = ceiling(nnoise/2)
nnC = nnoise-nnN
nd = genColumns(nrows,'n',nnN,nnC)
d = cbind(d,nd)
}
d$y = d$y > 0
d
}
# ------
# get the signal scores for the variables
# in a data set
# assume output is a binary variable named y
get_chiscores = function(dframe, varnames) {
nvar = length(varnames)
scores = numeric(nvar)
for(i in seq_len(nvar)) {
model = glm(paste("y~",varnames[i]), dframe,
family=binomial(link="logit"))
scores[i] = get_significance(model)
}
sframe = data.frame(var=varnames,
scores=scores, stringsAsFactors=FALSE)
sframe
}
#
# Plot the scores of each variable
# frm has columns var and scores
# (output of get_chiscores)
#
scoreplot = function(frm, threshold, sort=1) {
n = dim(frm)[1]
frm$var = reorder(frm$var, frm$scores*sort, FUN=sum)
frm$goodvar = frm$scores < threshold
ggplot(frm, aes(x=var, y=scores, ymin=0, ymax=scores, color=goodvar)) +
geom_pointrange() +
geom_hline(yintercept=threshold, color="red", linetype=2) +
scale_color_manual(values=c("TRUE"="darkgreen", "FALSE"="darkgray")) +
theme(legend.position="none")
}Example 1 Small example of evaluating a variable with signal, and without
set.seed(3266)
N = 1000
s1 = rnorm(N)
n1 = rnorm(N)
y = 2*s1 + rnorm(N)
dframe = data.frame(y=y>0, s1=s1, n1=n1)
nperm=500
score_variable(dframe, "y", "s1", nperm,
title='Signal variable deviance,')##
## Call:
## glm(formula = "y~x", family = binomial(link = "logit"), data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.80896 -0.42049 0.01506 0.45285 2.87230
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.02815 0.09930 -0.283 0.777
## x 3.50570 0.22463 15.606 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1386.23 on 999 degrees of freedom
## Residual deviance: 643.95 on 998 degrees of freedom
## AIC: 647.95
##
## Number of Fisher Scoring iterations: 6
## [1] "Chi-sq estimate: 1.91151650110488e-163"score_variable(dframe, "y", "n1", nperm,
title='Noise variable deviance,')##
## Call:
## glm(formula = "y~x", family = binomial(link = "logit"), data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.256 -1.183 1.116 1.170 1.243
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.01460 0.06329 0.231 0.818
## x -0.05703 0.06575 -0.867 0.386
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1386.2 on 999 degrees of freedom
## Residual deviance: 1385.5 on 998 degrees of freedom
## AIC: 1389.5
##
## Number of Fisher Scoring iterations: 3
## [1] "Chi-sq estimate: 0.385532666108562"Example 2
Small example of using chi-sq significance for variable filtering
ngood = 5; nnoise = 5; N = 1000
coefs = mkCoefs(ngood)
dframe = mkData(N, coefs, nnoise)
summary(dframe)## y gn_1 gn_2 gn_3
## Mode :logical Min. :-3.44603 Min. :-2.85755 Min. :-3.082419
## FALSE:648 1st Qu.:-0.67585 1st Qu.:-0.58603 1st Qu.:-0.662868
## TRUE :352 Median : 0.01202 Median : 0.02734 Median : 0.000484
## NA's :0 Mean :-0.01574 Mean : 0.04027 Mean :-0.015867
## 3rd Qu.: 0.63183 3rd Qu.: 0.74720 3rd Qu.: 0.678202
## Max. : 3.14300 Max. : 3.21589 Max. : 3.461750
## gc_1 gc_2 nn_1
## Length:1000 Length:1000 Min. :-3.15092
## Class :character Class :character 1st Qu.:-0.68175
## Mode :character Mode :character Median : 0.04855
## Mean :-0.03193
## 3rd Qu.: 0.61428
## Max. : 3.58665
## nn_2 nn_3 nc_1
## Min. :-2.94448 Min. :-2.87404 Length:1000
## 1st Qu.:-0.66348 1st Qu.:-0.73475 Class :character
## Median :-0.01824 Median :-0.05155 Mode :character
## Mean :-0.01459 Mean :-0.05146
## 3rd Qu.: 0.63712 3rd Qu.: 0.66751
## Max. : 3.84612 Max. : 3.26796
## nc_2
## Length:1000
## Class :character
## Mode :character
##
##
## varnames = setdiff(colnames(dframe), "y")
sframe = get_chiscores(dframe, varnames)
threshold = 0.05 # be generous in accepting a variable (1/20 false positive rate)
scoreplot(sframe, threshold)
print("Variables selected:")## [1] "Variables selected:"print(sframe[sframe$scores<threshold,])## var scores
## 2 gn_2 3.909215e-99
## 3 gn_3 3.561903e-09
## 4 gc_1 3.906063e-06
## 5 gc_2 4.325259e-06
## 8 nn_3 1.459305e-02print("True coefficients of signal variables")## [1] "True coefficients of signal variables"print(coefs)## $gn_1
## [1] 0.03406483
##
## $gn_2
## [1] -1.457936
##
## $gn_3
## [1] 0.4138757
##
## $gc_1
## $gc_1$a
## [1] -0.504306
##
## $gc_1$b
## [1] -0.2303952
##
## $gc_1$c
## [1] -1.174718
##
##
## $gc_2
## $gc_2$a
## [1] -0.2770224
##
## $gc_2$b
## [1] -0.1947211
##
## $gc_2$c
## [1] 0.2442891Example 3
Larger example; lots of noise
ngood = 5; nnoise = 2000; N = 2500
coefs = mkCoefs(ngood)
dframe = mkData(N, coefs, nnoise)
varnames = setdiff(colnames(dframe), "y")
print("True coefficients of signal variables")## [1] "True coefficients of signal variables"print(coefs)## $gn_1
## [1] 0.8907755
##
## $gn_2
## [1] 0.09630513
##
## $gn_3
## [1] 1.072262
##
## $gc_1
## $gc_1$a
## [1] 1.371897
##
## $gc_1$b
## [1] -0.8606408
##
## $gc_1$c
## [1] -0.4705882
##
##
## $gc_2
## $gc_2$a
## [1] -2.172888
##
## $gc_2$b
## [1] -0.9766713
##
## $gc_2$c
## [1] 0.9012387sframe = get_chiscores(dframe, varnames)
# we should always pick the threshold first, but
# we will use this as an example of what different
# thresholds will do
# 1/100 error rate, 1/40 error rate, 1/20 error rate
thresholds = c(0.01, 0.025, 0.05)
for(threshold in thresholds) {
n = sum(sframe$scores<threshold)
print(paste(n, "variables selected, threshold =", threshold))
print(sframe[sframe$scores<threshold,])
}## [1] "17 variables selected, threshold = 0.01"
## var scores
## 1 gn_1 4.744091e-40
## 3 gn_3 4.426174e-72
## 4 gc_1 3.754517e-64
## 5 gc_2 2.838399e-117
## 95 nn_90 7.278843e-03
## 353 nn_348 2.702906e-03
## 470 nn_465 3.625489e-03
## 617 nn_612 8.157911e-03
## 833 nn_828 7.904690e-03
## 1006 nc_1 2.983575e-03
## 1265 nc_260 6.512229e-03
## 1290 nc_285 9.559572e-03
## 1370 nc_365 5.091275e-03
## 1490 nc_485 7.549149e-03
## 1606 nc_601 2.689314e-03
## 1650 nc_645 5.622456e-03
## 1912 nc_907 9.193532e-03
## [1] "46 variables selected, threshold = 0.025"
## var scores
## 1 gn_1 4.744091e-40
## 3 gn_3 4.426174e-72
## 4 gc_1 3.754517e-64
## 5 gc_2 2.838399e-117
## 25 nn_20 1.534896e-02
## 50 nn_45 1.054503e-02
## 95 nn_90 7.278843e-03
## 98 nn_93 1.748200e-02
## 153 nn_148 1.463944e-02
## 338 nn_333 1.652431e-02
## 353 nn_348 2.702906e-03
## 386 nn_381 1.125829e-02
## 426 nn_421 1.037611e-02
## 470 nn_465 3.625489e-03
## 478 nn_473 1.465016e-02
## 503 nn_498 1.977930e-02
## 529 nn_524 2.162455e-02
## 563 nn_558 1.895795e-02
## 596 nn_591 1.226371e-02
## 617 nn_612 8.157911e-03
## 631 nn_626 2.439658e-02
## 665 nn_660 1.843377e-02
## 781 nn_776 2.499622e-02
## 833 nn_828 7.904690e-03
## 952 nn_947 2.143719e-02
## 978 nn_973 1.328706e-02
## 1006 nc_1 2.983575e-03
## 1059 nc_54 1.426079e-02
## 1187 nc_182 2.355118e-02
## 1230 nc_225 1.576357e-02
## 1265 nc_260 6.512229e-03
## 1290 nc_285 9.559572e-03
## 1370 nc_365 5.091275e-03
## 1490 nc_485 7.549149e-03
## 1522 nc_517 1.481355e-02
## 1549 nc_544 2.044036e-02
## 1606 nc_601 2.689314e-03
## 1615 nc_610 2.097747e-02
## 1650 nc_645 5.622456e-03
## 1810 nc_805 1.363454e-02
## 1840 nc_835 1.659934e-02
## 1847 nc_842 1.051913e-02
## 1893 nc_888 1.580964e-02
## 1912 nc_907 9.193532e-03
## 1918 nc_913 1.375041e-02
## 1998 nc_993 1.545178e-02
## [1] "103 variables selected, threshold = 0.05"
## var scores
## 1 gn_1 4.744091e-40
## 2 gn_2 3.591492e-02
## 3 gn_3 4.426174e-72
## 4 gc_1 3.754517e-64
## 5 gc_2 2.838399e-117
## 21 nn_16 3.219597e-02
## 25 nn_20 1.534896e-02
## 40 nn_35 3.404971e-02
## 42 nn_37 4.844582e-02
## 50 nn_45 1.054503e-02
## 53 nn_48 2.618896e-02
## 93 nn_88 3.240945e-02
## 95 nn_90 7.278843e-03
## 98 nn_93 1.748200e-02
## 144 nn_139 4.228818e-02
## 153 nn_148 1.463944e-02
## 244 nn_239 4.254657e-02
## 249 nn_244 4.379877e-02
## 282 nn_277 3.771266e-02
## 338 nn_333 1.652431e-02
## 353 nn_348 2.702906e-03
## 371 nn_366 4.976024e-02
## 383 nn_378 3.239405e-02
## 386 nn_381 1.125829e-02
## 388 nn_383 2.694754e-02
## 424 nn_419 3.732872e-02
## 426 nn_421 1.037611e-02
## 428 nn_423 2.653770e-02
## 433 nn_428 3.461036e-02
## 445 nn_440 2.518610e-02
## 460 nn_455 3.921065e-02
## 470 nn_465 3.625489e-03
## 478 nn_473 1.465016e-02
## 503 nn_498 1.977930e-02
## 524 nn_519 3.026417e-02
## 529 nn_524 2.162455e-02
## 552 nn_547 2.778526e-02
## 563 nn_558 1.895795e-02
## 565 nn_560 3.183238e-02
## 578 nn_573 4.064631e-02
## 596 nn_591 1.226371e-02
## 617 nn_612 8.157911e-03
## 631 nn_626 2.439658e-02
## 665 nn_660 1.843377e-02
## 697 nn_692 3.818386e-02
## 724 nn_719 3.046587e-02
## 726 nn_721 3.717127e-02
## 781 nn_776 2.499622e-02
## 784 nn_779 4.837252e-02
## 785 nn_780 2.672533e-02
## 832 nn_827 4.028992e-02
## 833 nn_828 7.904690e-03
## 848 nn_843 3.172189e-02
## 851 nn_846 3.228479e-02
## 852 nn_847 3.043253e-02
## 872 nn_867 3.829050e-02
## 885 nn_880 3.935873e-02
## 912 nn_907 2.991736e-02
## 952 nn_947 2.143719e-02
## 970 nn_965 3.839762e-02
## 978 nn_973 1.328706e-02
## 985 nn_980 2.555435e-02
## 1006 nc_1 2.983575e-03
## 1007 nc_2 3.160839e-02
## 1059 nc_54 1.426079e-02
## 1060 nc_55 3.073947e-02
## 1128 nc_123 2.984885e-02
## 1187 nc_182 2.355118e-02
## 1199 nc_194 2.514109e-02
## 1218 nc_213 2.816941e-02
## 1230 nc_225 1.576357e-02
## 1250 nc_245 4.320626e-02
## 1253 nc_248 4.925703e-02
## 1265 nc_260 6.512229e-03
## 1267 nc_262 2.641738e-02
## 1290 nc_285 9.559572e-03
## 1307 nc_302 2.732476e-02
## 1330 nc_325 4.497941e-02
## 1336 nc_331 3.163276e-02
## 1370 nc_365 5.091275e-03
## 1490 nc_485 7.549149e-03
## 1522 nc_517 1.481355e-02
## 1525 nc_520 4.977668e-02
## 1549 nc_544 2.044036e-02
## 1598 nc_593 2.518194e-02
## 1606 nc_601 2.689314e-03
## 1615 nc_610 2.097747e-02
## 1632 nc_627 2.792739e-02
## 1650 nc_645 5.622456e-03
## 1672 nc_667 3.779772e-02
## 1686 nc_681 3.638010e-02
## 1697 nc_692 4.903456e-02
## 1719 nc_714 3.363250e-02
## 1810 nc_805 1.363454e-02
## 1826 nc_821 2.870117e-02
## 1840 nc_835 1.659934e-02
## 1847 nc_842 1.051913e-02
## 1893 nc_888 1.580964e-02
## 1912 nc_907 9.193532e-03
## 1918 nc_913 1.375041e-02
## 1962 nc_957 3.220855e-02
## 1975 nc_970 2.818335e-02
## 1998 nc_993 1.545178e-02threshold = 0.01
scoreplot(sframe, threshold, sort=-1) + coord_flip()
Not try fitting a random forest on this data.
library('randomForest')## randomForest 4.6-10
## Type rfNews() to see new features/changes/bug fixes.# devtools::install_github("WinVector/WCPlots")
library('WVPlots')## Loading required package: grid
## Loading required package: gridExtra
## Loading required package: reshape2
## Loading required package: ROCR
## Loading required package: gplots
##
## Attaching package: 'gplots'
##
## The following object is masked from 'package:stats':
##
## lowess
##
## Loading required package: plyr
## Loading required package: stringr# convert all character vectors to factor for randomForest
dframe = as.data.frame(lapply(dframe, function(col) {
if(is.character(col)) {as.factor(col)}
else {col}
}))
allvars = sframe$var
selvars = sframe$var[sframe$scores < 0.01]
predframeTrain = data.frame(y=dframe$y)
# model on all variables
rfmodel_full = randomForest(x=dframe[,allvars], y=as.factor(dframe$y))
predframeTrain$trainfull = predict(rfmodel_full, newdata=dframe[,allvars], type='prob')[,'TRUE']
rfmodel_red = randomForest(dframe[,selvars], as.factor(dframe$y))
predframeTrain$trainred = predict(rfmodel_red, newdata=dframe[,selvars], type='prob')[,'TRUE']
testframe = mkData(N, coefs, nnoise)
testframe = as.data.frame(lapply(testframe, function(col) {
if(is.character(col)) {as.factor(col)}
else {col}
}))
predframeTest = data.frame(y=testframe$y)
predframeTest$testfull = predict(rfmodel_full, newdata=testframe[,allvars], type='prob')[,'TRUE']
predframeTest$testred = predict(rfmodel_red, newdata=testframe[,selvars], type='prob')[,'TRUE']
printutil = function(frame, column, yname, title) {
print(DoubleDensityPlot(frame, column, yname, title=title))
print(ROCPlot(frame, column, yname, title=title))
}
print("performance, all variables")## [1] "performance, all variables"printutil(predframeTrain, "trainfull", "y", "training performance, all variables")
printutil(predframeTest, "testfull", "y", "test performance, all variables")
print("performance, reduced variables")## [1] "performance, reduced variables"printutil(predframeTrain, "trainred", "y", "training performance, reduced variables")
printutil(predframeTest, "testred", "y", "test performance, reduced variables")
