Вопросы по R (статистика). Простой и посложнее.

R я не умею, но в матлабе у меня получилось 0.5448.

Так что правильно, по идее.

Спасибо, уже хоть что-то :)

когда хотят использовать результат фитинга модели говорят

predict(lm(data), newdata)

генерик predict() сам найдет метод

зачем «обезьяним», тут всё вручнуюТМ :)

> library(MASS)
> dose.p
function (obj, cf = 1:2, p = 0.5) 
{
    eta <- family(obj)$linkfun(p)
    b <- coef(obj)[cf]
    x.p <- (eta - b[1L])/b[2L]
    names(x.p) <- paste("p = ", format(p), ":", sep = "")
    pd <- -cbind(1, x.p)/b[2L]
    SE <- sqrt(((pd %*% vcov(obj)[cf, cf]) * pd) %*% c(1, 1))
    res <- structure(x.p, SE = SE, p = p)
    class(res) <- "glm.dose"
    res
}
<bytecode: 0x104ef78>
<environment: namespace:MASS>

Чего попросили, того и подогнала :)

> read.table("dat.txt", header=T)
      C   IC
1 1e+00 0.60
2 1e-01 0.30
3 1e-02 0.20
4 1e-03 0.03
5 1e-04 0.06
> data<-read.table("dat.txt", header=T)
> library(MASS)

> glm(IC~C,data=data,family=quasibinomial(link=probit))

Call:  glm(formula = IC ~ C, family = quasibinomial(link = probit), 
    data = data)

Coefficients:
(Intercept)            C  
     -1.108        1.401  

Degrees of Freedom: 4 Total (i.e. Null);  3 Residual
Null Deviance:	    1.196 
Residual Deviance: 0.3282 	AIC: NA



> plot(seq(from=0,to=1,length.out=25),  predict(glm(IC~C,data=data,family=quasibinomial(link=probit)), list(C=seq(from=0,to=1,length.out=25)),  type = "response"))

> points(seq(from=0,to=1,length.out=25),  predict(glm(IC~C,data=data,family=quasibinomial(link=probit)), list(C=seq(from=0,to=1,length.out=25)),  type = "response"), col="red")

Вот так оно подгоняет, вроде как синтетика сошлась :)

> data1 <- data.frame(doza=data$IC, dead <- data$C)
> plot(data1)
> points(seq(from=0,to=0.6,length.out=25),  predict(glm(dead~doza,data=data1,family=binomial(link="logit")), list(doza=seq(from=0,to=0.6,length.out=25)),  type = "response"), col="red")
Предупреждение
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!

> points(seq(from=0,to=0.6,length.out=25),  predict(glm(dead~doza,data=data1,family=binomial(link="probit")), list(doza=seq(from=0,to=0.6,length.out=25)),  type = "response"), col="green")
Предупреждение
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!


> dose.p(glm(dead~doza,data=data1,family=binomial(link="logit")))
              Dose        SE
p = 0.5: 0.3758331 0.2084593
Предупреждение
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
> 

> dose.p(glm(dead~doza,data=data1,family=binomial(link="probit")))              Dose        SE
p = 0.5: 0.3912484 0.1633893
Предупреждение
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
> 

Никак не пойму математический смысл коэффициента cf = 1:2.

Это номера коэффициентов модели имеющих отношение именно к дозе (поскольку модель может быть намного сложнее, например еще и пол и фазу луны учитывать)

вот там дословно так происходит, при таком значении

coef(model)[1:2]

Ясно, спасибо.

И вопрос ещё чисто математический. Сейчас только допёр, что пробит-распределение в целом виде уже содержит в себе и экспоненциальную, и линейную, и логарифмическую часть. Т.е. набор точек аппроксимируется фактически каким-то отрезком кривой такого распределения. Тогда имеет ли смысл брать логарифм от аргумента (в нашем случае C) с последующим экспоненцированием при вычислении IC50? Но с другой стороны результат, полученный таким образом выглядит на вид адекватнее, чем при пробит-регрессии без лишних преобразований...

Вообще, как всё запущенно у биологов :) В опыте по 4..6 измерений — уже много. Результаты как только могут не интерпретироваться...

1) В glm в виде link функции уже заложено преобразование данных (как правило линеаризующее зависимость).

Если осуществляется несколько последовательных преобразований, то лучше когда модель «знает о них», эта априорная информация часто позволяет производить подгонку оптимальным образом. Функция ошибок должна остаться одной и заключать в себя все особенности модели и наши априорные знания о ней.

2) Построенный по исходным синтетическим данным график показывает полную неадекватность модели.

Построенный график по исправленным синтетическим данным показывает что в нем просто нет никакой информации по заданному интервалу.

3) Сами синтетические данные не представляют из себя никакой «статистики», поскольку ни один опыт не был проведен множества раз. Если точка представляет из себя «свертку» группы животных получавших дозу, то отсутствия параметров группы не позволяют провести бутстреп оценку данной точки (или по всем группам бутстреп оценку для регрессионной модели в целом).

4) Результаты никак особенно интерпретироваться не могут, «окно интерпретаций» всегда крайне узкое. Такое впечатление образуется после курса теовера на физтехе, и оно в корне не верно (вообще я бы этот курс отменил или вынес на факультатив у технарей, вреда больше чем пользы).

Построенный по исходным синтетическим данным график показывает полную неадекватность модели.

Да, я примерно так и понял. Он лучше всего ложится на y=a*log(x)+b.

Кстати, как для glm посмотреть коэффициент корреляции?

Сами синтетические данные не представляют из себя никакой «статистики», поскольку ни один опыт не был проведен множества раз.

Это главная беда современной биохимии (не только российской, мировой :D). Меня, как классического химика их объёмы экспериментов просто убивают. Это нас учили заранее планировать количество серий, число экспериментов в серии... Тут же песец. В одном опыте как правило одна серия, число измерений единичное.

Если точка представляет из себя «свертку» группы животных получавших дозу

В используемых данных IC — это оценка ингибирования клеток, доля выживших через заданное время после воздействия препарата.

Такое впечатление образуется после курса теовера на физтехе, и оно в корне не верно

Ну, фиг знает, у нас был отдельный спецкурс «планирование эксперимента», кроме собственно курсов теорвера и курсов вычислительной математики :) Но методы, которые используются в биологии у нас и близко не упоминались.

На ту же тему: http://www.balancer.ru/g/p2256817

Сейчас только допёр, что пробит-распределение в целом виде уже содержит в себе и экспоненциальную, и линейную, и логарифмическую часть. Т.е. набор точек аппроксимируется фактически каким-то отрезком кривой такого распределения.

Там все эти экспоненты и логарифмы для другого используются. Пробит - это модель бинарного выбора, да/нет, с некоторой вероятностью. Если использовать обычную линейную регрессию, то она работать на такой задаче будет очень плохо. Поэтому используют специальную функцию, которая одним хвостом ассимптотически уходит в ноль, другим в единицу. А подготовка данных - отдельная задача. Я тоже сначала попробовал подогнать без логарифмирования, но результат был явно неадекватный, поскольку исходное распределение несимметрично. С логарифмированием - более чем.

1) посчитать коэффициент детерминации

https://stat.ethz.ch/pipermail/r-help/2010-June/243113.html

2) Вас там убивают параметрической статистикой :(

Никакая «нормальность» или соответствие другому «закону распределения» не нужны. Если есть экспериментальные данные, то у них есть дисперсия.

Все измерения вероятности (для любой сложности проверяемой гипотезы-модели) имея компьютер можно провести вычисляя напрямую по этой «экспериментальной дисперсии» с помощью вариантов МонтеКарло (рандомизация и перевыборки).

С помощью «законов распределения» и «критериев» считают что бы снизить объем вычислений необходимый для достижения нужной точности. Ну и что бы не городить сложные схемы счета для всяких крайних случаев.