Объяснение разницы в исходном и начальном значении

У меня есть следующая функция; (1) вычислить разницу отклонений для каждой переменной, которая у меня есть, и (2) выполнить начальную загрузку для разницы отклонений, которую я рассчитал на первом шаге.

set.seed(1001)

xfunction <- function(d,i)
{
glm.snp1 <- glm(PHENOTYPE~d[i], family="binomial", data=training1)
null <- glm.snp1$null.deviance
residual <- glm.snp1$deviance
dDeviance <- null-residual
return(dDeviance)
}

myboot <- function(d)
{
boot(d,xfunction, R=1000)
}

result <- lapply(training1,function(x)myboot(x))

Таким образом, в основном из результата я получил значения для исходного dDeviance (без начальной загрузки), и я могу рассчитать среднее значение (dDeviance) из начальной загрузки. Мне нужна помощь в объяснении, почему исходное и среднее значение бутстрапа слишком разные? Например, для одной из переменных исходное значение dDeviance равно 0,024, а среднее значение начальной загрузки dDeviance равно 0,000412.


r statistics-bootstrap resampling
person Shima    schedule 24.02.2016    source источник
comment
Это кажется маловероятным, поскольку распределение бутстрапированной статистики в основном соответствует статистике в исходных данных. Глядя на ваш код, я думаю, что ваша проблема заключается в индексации. Вы используете d[i], но по-прежнему используете исходный обучающий набор и исходный фенотип.   -  person Heroka    schedule 24.02.2016
comment
Спасибо @Heroka за комментарий. Я тоже согласен с вами, что значение должно быть около исходного значения. Я пытался избавиться от данных в xfunction, но он не может определить, откуда взялся мой ФЕНОТИП. Я признателен, если вы могли бы предложить поправки в моей функции. К вашему сведению, мой результат - ФЕНОТИП, а предикторы состоят из 1500 SNP с 1000 образцами.   -  person Shima    schedule 24.02.2016
comment
В моем опыте с начальной загрузкой у меня была функция, которая брала набор данных, содержащий все, что мне нужно, и индексы, а затем использовала это. Таким образом, фенотип должен быть на тренировке1, тогда вы сможете подогнать свою модель под data=d[i,].   -  person Heroka    schedule 24.02.2016
comment
Я согласен с @Heroka, вы должны подмножить data.frame во время начальной загрузки, а не в формуле. Я тоже не понимаю, что ты там делаешь со своим циклом lapply.   -  person Roland    schedule 24.02.2016
comment
Спасибо вам обоим за комментарии. Я постараюсь над этим поработать. Ну, идея, которую я поместил туда, заключается в том, что он вычисляет dDeviance для каждого столбца в training1[,2:1501] и загружает значение с помощью функции загрузки.   -  person Shima    schedule 24.02.2016

Ответы (1)


arrow_upward
0
arrow_downward

Как указано в комментариях, лучше подмножить индексы data.frame, чтобы сделать их совместимыми с boot. Если вы хотите выполнить итерацию по различным переменным и применить эту функцию, лучше сделать это внутри итерации.

Для каждого бутстрапа вы подбираете разные модели и извлекаете значения.

Итак, мы переписываем функцию, и она принимает дополнительный параметр формулы, indep_var, который указывает столбцы для регрессии PHENOTYPE:

xfunction <- function(d,ind,indep_var){

    dDeviance = sapply(indep_var,function(V){
         f = reformulate(V,response="PHENOTYPE")  
         glm.snp1 <- glm(f, family="binomial", data=d[ind,])
         glm.snp1$null.deviance-glm.snp1$deviance
    })
return(dDeviance)
}

Мы можем настроить пример набора данных:

set.seed(111)
training1 = data.frame(PHENOTYPE=rbinom(100,1,0.5),matrix(rnorm(500),ncol=5))

head(training1)
  PHENOTYPE         X1          X2         X3         X4         X5
1         1  0.1916634 -0.09152026 -0.9685094 -1.0503824 -0.9137690
2         1  1.5525443 -1.87430581  0.9119331  0.3251424  0.1126909
3         0  0.9142423 -0.66416202  0.0928326 -2.1048716 -2.3597249
4         1  0.3586254  0.20341282 -0.5329309 -0.9551027 -1.5693983
5         0  0.1750956 -2.59444339 -1.6669055 -0.5306399  1.2274569
6         0 -0.8472678 -0.09373393 -0.5743455  0.8274405  0.7620480

И попробуйте это на данных:

xfunction(training1,1:nrow(training1),indep_var=c("X1","X2","X3"))
       X1        X2        X3 
0.1189847 0.2512539 0.2684927

Затем с помощью бутстрапа:

library(boot)
boot(training1,xfunction,R=1000,indep_var=c("X1","X2","X3"))

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = training1, statistic = xfunction, R = 1000, indep_var = c("X1", 
    "X2", "X3"))


Bootstrap Statistics :
     original   bias    std. error
t1* 0.1189847 1.033286    1.564971
t2* 0.2512539 1.047503    1.863012
t3* 0.2684927 1.005062    1.747959
person StupidWolf    schedule 30.08.2020