Kodomo

Пользователь

Учебная страница курса биоинформатики,
год поступления 2014

Линейные модели часть2

Презентация тут: [лекция 7]

Сравнение моделей с помощью ANOVA

laptop=read.csv(file="laptop_price.csv", header=T)
fit2=lm(Price_RUR ~ Memory_Gb+HDD_Gb+HDD_type, data=laptop)
fit1=lm(Price_RUR ~ Memory_Gb, data=laptop)
anova(fit1, fit2)

Противоположный пример 

fit1=lm(Price_RUR ~ Memory_Gb+HDD_Gb+HDD_type+Color,
          data=laptop)
fit2=lm(Price_RUR ~ Memory_Gb+HDD_Gb+HDD_type, data=laptop)
anova(fit1, fit2)

Предсказание с помощью модели

L3 = lm(formula = Price_RUR ~ HDD_Gb + HDD_type + HDD_Gb:HDD_type, data = laptop)
newlaptops=data.frame( HDD_Gb=c(200, 1000,                               500), HDD_type=c("SSD", "HDD", "HDD"))
predict( L3, newlaptops )

Кросс-валидация

library(DAAG)
cv.lm(laptop, L3, m=5)

Дополнительный аргумент – функция cost, по умолчанию: 

cost=function(y, ypred) { mean ( (y-ypred)^2 ) }

glm – обобщенные линейные модели

Регрессия Пуассона

Пример: растет ли посещаемость сайта со временем gaData.rda 

load('gaData.rda')
gaData$julian<-julian(gaData$date)

plot(gaData$visits~gaData$julian, xlab="julian",      ylab="visits", pch=19, col="darkgrey")
lm1<-lm(gaData$visits~gaData$julian)
abline(lm1, col='red', lwd=3)
glm1<-glm(gaData$visits~gaData$julian, family='poisson')
lines(gaData$julian, glm1$fitted, col='blue', lwd=3)

Логистическая регрессия

Пример: таблица выигрышей команды. Хотим вычислять вероятность выигрыша в зависимости от количества очков ravensData.rda 

load("ravensData.rda")

lmRav<-  lm(ravensData$ravenWinNum~ravensData$ravenScore)
plot(ravensData$ravenScore, lmRav$fitted, pch=19, col='red', ylab="Prob win")

logReg<-glm(  ravensData$ravenWinNum~ravensData$ravenScore,   family="binomial")
logReg

plot(ravensData$ravenScore, ravensData$ravenWinNum,      pch=19, col='blue')

points(ravensData$ravenScore, logReg$fitted, pch=19,        col='red')

Кросс-валидация для категориальной зависимой переменной

Предскажем наличие SSD в ноутбуке по его цене и объему диска 

gl1=glm(HDD_type ~ Price_RUR + HDD_Gb,
          family=binomial, data=laptop)
gl1
plot(predict( gl1 ) ~ laptop$HDD_type)
plot(predict( gl1, type="response" ) ~ laptop$HDD_type)

cost=function( y, ypred ) {
  mean( abs(y-ypred)>0.5 )
}

library(boot)
cv1=cv.glm(laptop, gl1, cost, K=3)

Задача

Исходные данные. в аннотированных энхансерах и промоторах в геноме вычислены уровни различных эпигенетических модификаций: модификации гистонов H3K27ac, H3K27me3, H3K36me3, H3K4me1, H3K4me3, H3K9me3, ChIP-seq на РНК-полимеразу II (Pol2) и CTCF.

В таблице http://makarich.fbb.msu.ru/artemov/R/enhancer_promoter.tab приведены уровни эпигенетических сигналов, а первый столбец Type указывает на то, является отрезок энхансером (1) или промотором (0).

Задача: научиться отличать по эпигенетическим модификациям промоторы от энхансеров. 

tab2=read.table("http://makarich.fbb.msu.ru/artemov/R/enhancer_promoter.tab", header=T, stringsAsFactors=F)

Визуализация данных

library(gplots)

#Корреляции между предикторами, а также между предикторами и outcome
heatmap.2(cor( tab2 ), col="bluered", symbreaks=T)

heatmap: возьмем 100 случайных промоторов и 100 случайных энхансеров: 

heatmap.2(data.matrix(tab2[c(sample(which(tab2$Type==0), 100), sample(which(tab2$Type==1), 100)),2:ncol(tab2)]), trace="none", col="bluered", scale="column", breaks= (-20:20)/10, Rowv=F, RowSideColors=c(rep("red", 100), rep("orange", 100)))

#более детально:
#вектор номеров строк, первые 100 указывают на строки с промоторами, последние 100 - на строки с энхансерами
#c(sample(which(tab2$Type==0), 100), sample(which(tab2$Type==1), 100))
#Не перемешиваем по строкам (уже упорядочены так, как нам надо): Rowv=F, 
#Цвет сбоку от heatmap-а указывает, является элемент промотором (красный) или энхансером (оранжевый):
#RowSideColors=c(rep("red", 100), rep("orange", 100))

Построение и сравнение обобщенных линейных моделей (логистическая регрессия)

Для начала, проведем логистическую регрессию от одной переменной, например, H3k4me3. По тем коэффициентам регрессии, которые получаются после обучения, модель строит линейную комбинацию предикторов - число от -Inf до Inf. Затем данное число нужно превратить в вероятность, что Type==1 (то есть объект - энхансер, а не промотер). Такие вероятности хранятся в векторе gl$fitted (для тех точек, на которых модель была обучена). 

gl=glm(Type ~ H3k4me3, family=binomial, data=tab2)
subsamp=sample(1:nrow(tab2), 200)
plot(tab2[subsamp,]$H3k4me3, tab2[subsamp,]$Type, pch=19, col='blue', xlab="H3K4me3", ylab="Probability of being an enhancer")
points(tab2[subsamp,]$H3k4me3, gl$fitted[subsamp], pch=19, col='red')
legend("topright", legend=c("Real (0 or 1)", "Predicted by model"), fill=c("blue", "red"))

Модель с двумя предикторами (H3k4me1 + H3k4me3) 

gl2=glm(Type ~ H3k4me1 + H3k4me3, family=binomial, data=tab2)
anova(gl2, gl, test="Chisq")

Сравнение моделей 

gl1=glm(Type ~ ., family=binomial, data=tab2)
gl2=glm(Type ~ H3k4me1 + H3k4me3, family=binomial, data=tab2)
anova(gl2, gl1)
anova(gl2, gl1, test="Chisq")

Кросс-валидация. 

cost=function( y, ypred ) {
  mean( abs(y-ypred)>0.5 )
}

library(boot)
cv1=cv.glm(tab2, gl, cost, K=5)
cv1$delta

Примените кросс-валидацию к другим возможным моделям: например, к моделям с одним предиктором или к модели, учитывающей все предикторы (Type ~ .).

Можно ли без ущерба оставить только два предиктора?

(*) Автоматический выбор параметров 

library(glmulti)
gl=glmulti(Type ~ ., family=binomial, data=tab2, level=1)

To be continued (SVM и другие классификаторы).