Master R 6 - Beyond the linear trend line

The modeling workflow
1 First, fit the model with the main predictors and all the relevant confounders, and then reduce the number of confounders by dropping out the non-significant ones. 
2 Decide whether to use the continuous variables in their original or categorized form.
3 Try to achieve a better fit by testing for non-linear relationships, if they are pragmaticaly relevant.
4 Finally check the model assumptions.

Logistic regression
library(catdata)
data(deathpenalty)
library(vcdExtra)
deathpenalty.expand <- expand.dft(deathpenalty)
binom.model.0 <- glm(DeathPenalty ~ DefendantRace, data = deathpenalty.expand, family = binomial)
summary(binom.model.0)
exp(cbind(OR = coef(binom.model.0), confint(binom.model.0)))

binom.model.1 <- update(binom.model.0, . ~ . + VictimRace)
summary(binom.model.1)
exp(cbind(OR = coef(binom.model.1), confint(binom.model.1)))

prop.table(table(factor(deathpenalty.expand$VictimRace,labels = c("VictimRace=0", "VictimRace=1")),factor(deathpenalty.expand$DefendantRace, labels = c("DefendantRace=0", "DefendantRace=1"))), 1)

Data consideration
A general rule of thumb, logistic regression models require at least 10 events per predictors.

Goodness of model fit
library(lmtest)
lrtest(binom.model.1)
library(BaylorEdPsych)
PseudoR2(binom.model.1)

Model comparison
lrtest(binom.model.0, binom.model.1)

Models for count data
Poisson regression
dfa <- readRDS('SMART_2013.RData')
(ct <- xtabs(~model+failure, data=dfa))
dfa <- dfa[dfa$model %in% names(which(rowSums(ct) - ct[, 1] > 0)),]

library(ggplot2)
ggplot(rbind(dfa, data.frame(model='All', dfa[, -1] )), aes(failure)) + ylab("log(count)") + 
geom_histogram(binwidth = 1, drop=TRUE, origin = -0.5)  + 
scale_y_log10() + scale_x_continuous(breaks=c(0:10)) + 
facet_wrap( ~ model, ncol = 3) + ggtitle("Histograms by manufacturer") + theme_bw()

poiss.base <- glm(failure ~ model, offset(log(freq)),  family = 'poisson', data = dfa)
summary(poiss.base)
contrasts(dfa$model, sparse = TRUE)
lrtest(poiss.base)

Negative binomial regression
library(MASS)
model.negbin.0 <- glm.nb(failure ~ model, offset(log(freq)), data = dfa)
lrtest(poiss.base,model.negbin.0)

Multivariate non-linear models
model.negbin.1 <- update(model.negbin.0, . ~ . + capacity_bytes + age_month + temperature)
model.negbin.2 <- update(model.negbin.1, . ~ . + PendingSector)
lrtest(model.negbin.0, model.negbin.1, model.negbin.2)
summary(model.negbin.2)
exp(data.frame(exp_coef = coef(model.negbin.2)))
dfa$model <- relevel(dfa$model, 'WDC')


model.negbin.3 <- update(model.negbin.2, data = dfa)
library(broom)
format(tidy(model.negbin.3), digits = 4)

library(data.table)
dfa <- data.table(dfa)
dfa[, temp6 := cut2(temperature, g = 6)]
temperature.weighted.mean <- dfa[, .(wfailure = weighted.mean(failure, freq)), by = temp6] 
ggplot(temperature.weighted.mean, aes(x = temp6, y = wfailure)) +  
geom_bar(stat = 'identity') + xlab('Categorized temperature') + ylab('Weighted mean of disk faults') + theme_bw()

model.negbin.4 <- update(model.negbin.0, .~. + capacity_bytes + age_month + temp6 + PendingSector, data = dfa)
AIC(model.negbin.3,model.negbin.4)

weighted.means <- rbind(dfa[, .(l = 'capacity', f = weighted.mean(failure, freq)), by = .(v = capacity_bytes)], dfa[, .(l = 'age', f = weighted.mean(failure, freq)), by = .(v = age_month)])
ggplot(weighted.means, aes(x = l, y = f)) + geom_step() +
facet_grid(. ~ v, scales = 'free_x') + theme_bw() +
ylab('Weighted mean of disk faults') + xlab('')

dfa[, capacity_bytes := as.factor(capacity_bytes)]
dfa[, age8 := cut2(age_month, g = 8)]
model.negbin.5 <- update(model.negbin.0, .~. + capacity_bytes + age8 + temp6 + PendingSector, data = dfa)
AIC(model.negbin.5, model.negbin.4)
format(tidy(model.negbin.5), digits = 3)

tmnb5 <- tidy(model.negbin.5)
str(terms <- tmnb5$term[tmnb5$p.value < 0.05][-1])

library(plyr)
ci <- ldply(terms, function(t) confint(model.negbin.5, t))
names(ci) <- c('min', 'max')
ci$term <- terms
ci$variable <- sub('[A-Z0-9\\]\\[,() ]*$', '', terms, perl = TRUE)

ggplot(ci, aes(x = factor(term), color = variable)) + 
geom_errorbar(ymin = min, ymax = max) + xlab('') +
ylab('Coefficients (95% conf.int)') + theme_bw() + 
theme(axis.text.x = element_text(angle = 90, hjust = 1), legend.position = 'top')
PseudoR2(model.negbin.6 )

Logistic Regression Model Calibration in R

all <- fread("all.txt", head=F, sep='\t')
names=c("HasClick", "PredScore", "Score");
setnames(all, names);
setkey(all, Score);
all[HasClick==-1, HasClick:=0];
all[,list(Impressions=.N, CTR=sum(HasClick)/.N), by=HasClick][order(-CTR)];
summary(all$HasClick)

#PredScore
all[PredScore==-1, PredScore:=0];
summary(all$PredScor)
pred<-prediction(all$PredScore,all$HasClick);
auc.tmp <- performance(pred,"auc");
auc <- as.numeric(auc.tmp@y.values);
auc
#[1] 0.6567984
#[1] /homes/lingh/all.txt
perf <- performance(pred,"tpr","fpr");
# KS is the maximum difference between the cumulative true positive and cumulative false positive rate.
max(attr(perf,'y.values')[[1]]-attr(perf,'x.values')[[1]]);
#[1] 0.244848

#p Score
all$Score <- all$Score
all$p <- 1/(1+exp(-all$Score))
summary(all$p)
pred<-prediction(all$p,all$HasClick);
auc.tmp <- performance(pred,"auc");
auc <- as.numeric(auc.tmp@y.values);
auc
#[1] 0.7015954
#[1] 0.6919405
#[1] 0.6992571
#[1] 0.7176237
perf <- performance(pred,"tpr","fpr");
# KS is the maximum difference between the cumulative true positive and cumulative false positive rate.
max(attr(perf,'y.values')[[1]]-attr(perf,'x.values')[[1]]);
#[1] 0.3135409
#[1] 0.3065308
#[1] 0.3129876
#[1] 0.3337304


# prior for bias correction
totimp <- nrow(all)
totclk <- sum(all$HasClick==1)
t <- totclk/totimp; t
p <-mean(all$p); p
all$p2 <- 1/(1+exp(-all$Score+log((1-t)*p/(t*(1-p)))));
mean(all$p2)
summary(all$p2)
pred<-prediction(all$p2,all$HasClick);
auc.tmp <- performance(pred,"auc");
auc <- as.numeric(auc.tmp@y.values);
auc
#[1] 0.7015954
#[1] 0.6919405
[1] 0.6992571
perf <- performance(pred,"tpr","fpr");
# KS is the maximum difference between the cumulative true positive and cumulative false positive rate.
max(attr(perf,'y.values')[[1]]-attr(perf,'x.values')[[1]]);
#[1] 0.3135409
#[1] 0.3065308
[1] 0.3129876


# ensemble
all$p3 <- (all$p2+all$PredScore)/2
summary(all$p3)
pred<-prediction(all$p3,all$HasClick);
auc.tmp <- performance(pred,"auc");
auc <- as.numeric(auc.tmp@y.values);
auc
#[1] 0.6996736
#[1] 0.6784338
[1] 0.6973308
perf <- performance(pred,"tpr","fpr");
# KS is the maximum difference between the cumulative true positive and cumulative false positive rate.
max(attr(perf,'y.values')[[1]]-attr(perf,'x.values')[[1]]);
#[1] 0.3098282
#[1] 0.2729899
[1] 0.3105527

clkinbin <- round(totclk/100); clkinbin
clkinbin <- 5
grp <- 1
tmp <- rep(1, totimp)
cnt <- 0
for (i in 1:totimp){
#for (i in 1:100000){
  if (all$HasClick[i]==1){
    cnt=cnt+1;
  }
  if (cnt==clkinbin){
    grp=grp+1;
    cnt=0;
  }
  tmp[i]=grp
}
all$grp=tmp
table(all$grp)
#all[grp==108, grp:=107];

groups=all[,list(Impressions=.N, Clicks=sum(HasClick), CTR=sum(HasClick)/.N,
                 AvgScore1=mean(PredScore),
                 AvgScore2=mean(p),
                 AvgScore3=mean(p2),
                 AvgScore4=mean(p3)),        
by=grp][order(grp)];
#write.csv(groups, "junk.csv")
with(groups, cor(CTR, AvgScore1))
#[1] 0.8853131
with(groups, cor(CTR, AvgScore2))
#[1] 0.87801
with(groups, cor(CTR, AvgScore3))
#[1] 0.8791485
with(groups, cor(CTR, AvgScore4))
# [1] 0.8803879

data <- groups
ggpairs(data[, c("CTR", "AvgScore1", "AvgScore2", "AvgScore3", "AvgScore4"), with=F])

sum((groups$AvgScore1-groups$CTR)^2)/length(groups$AvgScore1)
#[1] 2.920109e-06
sum((groups$AvgScore2-groups$CTR)^2)/length(groups$AvgScore2)
#[1] 0.0001156954
sum((groups$AvgScore3-groups$CTR)^2)/length(groups$AvgScore3)
#[1] 0.0001410945
#[1] 0.0001240449
sum((groups$AvgScore4-groups$CTR)^2)/length(groups$AvgScore4)
#[1] 4.574389e-05
#[1] 4.09054e-05


data<-groups
ggplot(data, aes(y=data$CTR, x=data$grp)) +
  geom_line(data=data, aes(y=data$CTR, x=data$grp), color='black', linetype=1, size = 1)  +
  geom_line(data=data, aes(y=data$AvgScore1, x=data$grp), color='red', linetype=1, size = 1)  +
#  geom_line(data=data, aes(y=data$AvgScore2, x=data$grp), color='purple', linetype=1, size = 1)  +
  geom_line(data=data, aes(y=data$AvgScore3, x=data$grp), color='darkgreen', linetype=1, size = 1)  +
  geom_line(data=data, aes(y=data$AvgScore4, x=data$grp), color='darkblue', linetype=1, size = 1)  +
  scale_x_continuous(limits=c(1, 100)) +
  scale_y_continuous(limits=c(0, 0.005)) +
  ggtitle("Binning Avg of CTR, HModel and VWModel") +
  xlab("Bin") + ylab("CTR") +
  theme(plot.title = element_text(face = "bold", size = 20)) +
  theme(axis.text.x = element_text(face = "bold", size = 16)) +
  theme(axis.text.y = element_text(face = "bold", size = 16)) +
  theme(axis.title.x = element_text(face = "bold", size = 16)) +
  theme(axis.title.y = element_text(face = "bold", size = 16, angle = 90)) +
  theme(legend.position = "top") +
  theme(legend.key = element_rect(colour = NA)) +
  theme(legend.title = element_blank())


# fit linear model
data<-groups[1:107]
with(data, cor(CTR, AvgScore1))
with(data, cor(CTR, AvgScore3))
m <- lm(CTR ~ AvgScore3, data = groups)
data$tmp <- 0.0005360169 + 0.0717973293*data$AvgScore3
sum((data$tmp-data$CTR)^2)/length(data$tmp)
# [1] 2.400827e-07

# fit tobit model
summary(m <- vglm(CTR ~ AvgScore3, tobit(Upper = 0.05), data = groups))
ctable <- coef(summary(m))
pvals <- 2 * pt(abs(ctable[, "z value"]), df.residual(m), lower.tail = FALSE)
cbind(ctable, pvals)
data$tmp <- 0.0005360131 + 0.0717974011*data$AvgScore3  
sum((data$tmp-data$CTR)^2)/length(data$tmp)


# fit loess smoother
m <- loess(CTR ~ AvgScore3, data = groups)
groups$tmp <- m$fitted
sum((groups$tmp-groups$CTR)^2)/length(groups$tmp)
# [1] 1.765797e-07


# fit Iso
install.packages("Iso")
require(Iso)
z = ufit(groups$CTR, groups$AvgScore3, type='b')
z = pava(groups$AvgScore3)
groups$tmp <- z
sum((groups$tmp-groups$CTR)^2)/length(groups$tmp)
#0.0001410945

data <- groups[1:105]
ggplot(data, aes(y=data$CTR, x=data$AvgScore3)) +
  geom_point(data=data, aes(y=data$CTR, x=data$AvgScore3), color='blue', shape=3, size = 3)  +
  geom_line(data=data, aes(y=data$AvgScore3, x=data$AvgScore3), color='black', linetype=3, size = 1)  +
  scale_x_continuous(limits=c(0, 0.0015)) +
  scale_y_continuous(limits=c(0, 0.0015)) +
  ggtitle("Distribution of CTR Vs LRScore") +
  xlab("LRScore") + ylab("CTR") +
  theme(plot.title = element_text(face = "bold", size = 20)) +
  theme(axis.text.x = element_text(face = "bold", size = 16)) +
  theme(axis.text.y = element_text(face = "bold", size = 16)) +
  theme(axis.title.x = element_text(face = "bold", size = 16)) +
  theme(axis.title.y = element_text(face = "bold", size = 16, angle = 90)) +
  theme(legend.position = "top") +
  theme(legend.key = element_rect(colour = NA)) +
  theme(legend.title = element_blank())

groups$nonClicks <-  groups$Impressions-groups$Clicks
zfit <- glm(cbind(groups$Clicks, groups$nonClicks ) ~ + data1$AvgScore3 , family = binomial, data=groups);
summary(lrfit);

glm.out = glm(cbind(Clicks, nonclick) ~ AvgScore2, family=binomial(logit), data=data)
plot(Clicks/Impressions ~ AvgScore2, data=data)
lines(data$AvgScore2, glm.out$fitted, type="l", col="red")
title(main="Group Data with Fitted Logistic Regression Line")

all$p3<-1/(1+exp(-all$Score*7.16219+8.82040))
summary(all$p3)
pred<-prediction(all$p3,all$HasClick);
auc.tmp <- performance(pred,"auc");
auc <- as.numeric(auc.tmp@y.values);
auc
#[1] 0.7015954
perf <- performance(pred,"tpr","fpr");
# KS is the maximum difference between the cumulative true positive and cumulative false positive rate.
max(attr(perf,'y.values')[[1]]-attr(perf,'x.values')[[1]]);
#[1] 0.3135409

Logistic Regression with Penalty Factor in R

The main goal of penalty is to minimize the sum of the unpenalized cost function plus the penalty term, which can be understood as adding bias and preferring a simpler model to reduce the variance in the absence of sufficient training data to fit the model.

- The ridge penalty shrinks the coefficients of correlated predictors towards each other (L2)
- The lasso tends to pick one of them and discard the others (L1)
- The elastic-net penalty mixes these two.

If predictors are correlated in groups, an $\alpha =0.5$ tends to select the groups in or out together. This is a higher level parameter, and users might pick 
a value upfront, else experiment with a few different values. 
One use of α is for numerical stability; for example, the elastic net with α=1−ϵ for some small ϵ>0 performs much like the lasso, but removes any degeneracies and wild behavior caused by extreme correlations.

## http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html#log
install.packages("glmnet", repos = "http://cran.us.r-project.org")
require(glmnet)

load("BinomialExample.RData")
fit = glmnet(x, y, family = "binomial")
summary(fit)
plot(fit, xvar = "dev", label = TRUE)

predict(fit, newx = x[1:5,], type = "class", s = c(0.05, 0.01))
cvfit = cv.glmnet(x, y, family = "binomial", type.measure = "class")
plot(cvfit)
cvfit$lambda.min
cvfit$lambda.1se
coef(cvfit, s = "lambda.min")

predict(cvfit, newx = x[1:10,], s = "lambda.min", type = "class")

Logistic Regression in R

Logistic regression relates a binary outcome such as purchase to predictors that may include continuous and factor variables, by modeling the variables' association with the probability of outcomes. 

A logistic regression model, also known as a logit model, is a member of the generalized linear models family. 

Coefficients in a logit model can be interpreted in terms of odds ratios, the degree to which they are associated with the increased or decrease likelihood of an outcome. 

A statistically significant result does not always mean that the model is appropriate. It is important to explore data thoroughly and to construct models on the basis of careful consideration. 

#########################################################
Fit with ratio and weights
#########################################################
data$network=as.factor(data$network)
data$Type=as.factor(data$Type)
model <- glm(Avg...Viewed ~Type, weights = Total.Hrs.Viewed, family=binomial, data=data)
summary(model)

## odds ratios only
exp(coef(model))
## CIs using profiled log-likelihood
confint(model)
## CIs using standard errors
confint.default(model)

wald.test(b = coef(model), Sigma = vcov(model), Terms = 2)

#########################################################
Fit with logical vertor of a two-level factor
#########################################################
hive -e 'SELECT freq1d, count(*), count(distinct user_id), sum(click), sum(action) from tmp5 group by freq1d order by freq1d;' > /homes/result1d.csv;

data <- read.table("result1d.csv", head =F);

names(data) <- c("freq", "imp", "uuser", "click", "action");

data$CTR <- data$click/data$imp;

data$CVR <- data$action/data$imp;

data$freq[data$action==6]

data1 <- data[1:49,];

data1$nonaction <- data1$imp-data1$action;

data1$freq1d <- as.factor(data1$freq);

lrfit <- glm(cbind(data1$action, data1$nonaction) ~ + data1$freq1d , family = binomial, data=data1);

summary(lrfit);

#########################################################

Fit downsampling data

#########################################################

train <- fread("train.csv", head=F, sep=',', colClasses=c("character", "character", "integer", "integer", "integer", "integer", "integer", "integer", "integer", "integer", "integer", "integer", "integer"));

test <- fread("test.csv", head=F, sep=',', colClasses=c("character", "character", "integer", "integer", "integer", "integer", "integer", "integer", "integer", "integer", "integer", "integer", "integer"));

setnames(train, c("rowid", "userid", "impdate", "imphour", "click", "action", "freq1d", "freq3d", "freq5d", "freq1w", "freq2w", "freq3w", "rect"));

setnames(test, c("rowid", "userid", "impdate", "imphour", "click", "action", "freq1d", "freq3d", "freq5d", "freq1w", "freq2w", "freq3w", "rect"));

train1 <- (train[ ,6:13, with=F]);

test1 <- (test[,6:13, with=F]);

train1$action <- factor(train1$action);

test1$action <- factor(test1$action);

# simple model;

lrfit <- glm(train1$action ~ ., family = binomial, data=train1);

summary(lrfit)

test1$score<-predict(lrfit,type='response',test1)

#calculate auc;

pred<-prediction(test1$score,test1$action);

auc.tmp <- performance(pred,"auc"); 

auc <- as.numeric(auc.tmp@y.values);

auc;

perf <- performance(pred,"tpr","fpr");

plot(perf);

# KS is the maximum difference between the cumulative true positive and cumulative false positive rate.

max(attr(perf,'y.values')[[1]]-attr(perf,'x.values')[[1]])

# top 3 variables

#works by sorting coefficient terms in equation and selecting top 3 in sort for each loan scored

g<-predict(lrfit,type='terms',test1)

ftopk<- function(x,top=3){

  res=names(x)[order(x, decreasing = TRUE)][1:top]

  paste(res,collapse=";",sep="")

}

topk=apply(g,1,ftopk,top=3)

test1<-cbind(test1, topk)

# add rect interaction;

lrfit <- glm(train1$action ~ freq1d+freq3d+freq5d+freq1w+freq2w+freq3w+rect+freq1d:rect+freq3d:rect+freq5d:rect+freq1w:rect+freq2w:rect+freq3w:rect, family = binomial, data=train1);

summary(lrfit)

test1$score<-predict(lrfit,type='response',test1)

#calculate auc;

pred<-prediction(test1$score,test1$action);

auc.tmp <- performance(pred,"auc"); 

auc <- as.numeric(auc.tmp@y.values); auc;

perf <- performance(pred,"tpr","fpr");

plot(perf);

# KS is the maximum difference between the cumulative true positive and cumulative false positive rate.

max(attr(perf,'y.values')[[1]]-attr(perf,'x.values')[[1]])