Использование весов в полиномиальной GAM (mgcv) - PullRequest
Новая
       17

Использование весов в полиномиальной GAM (mgcv)

0 голосов
/ 09 мая 2018

Мне интересно понять факторы, которые влияют на долю рыбы, которая покидает реку в виде мальков, парра или смолта (сумма равна 1). Я знаю, что это можно сделать с помощью nnet :: multinom, который прекрасно работает, но я хотел посмотреть, смогу ли я согласовать эти отношения со сплайнами, используя пакет mgcv. Чтобы сделать это, я собирался подогнать полиномиальную модель с 3 классами с весами регрессии относительно доли этого класса, наблюдаемой в тот день. Я могу приспособить полиномиальную модель с mgcv, однако, она не принимает веса. Ниже приведен пример - выходные данные идентичны с весами и без них. Кто-нибудь знает, возможно ли использовать весовые коэффициенты в полиномиальной игре, реализованной в mgcv? 

library(mgcv)
set.seed(6)
## simulate some data from a three class model
n <- 1000
f1 <- function(x) sin(3*pi*x)*exp(-x)
f2 <- function(x) x^3
f3 <- function(x) .5*exp(-x^2)-.2
f4 <- function(x) 1
x1 <- runif(n);x2 <- runif(n)
eta1 <- 2*(f1(x1) + f2(x2))-.5
eta2 <- 2*(f3(x1) + f4(x2))-1
p <- exp(cbind(0,eta1,eta2))
p <- p/rowSums(p) ## prob. of each category 
cp <- t(apply(p,1,cumsum)) ## cumulative prob.
## simulate multinomial response with these probabilities
## see also ?rmultinom
y <- apply(cp,1,function(x) min(which(x>runif(1))))-1
## plot simulated data...
plot(x1,x2,col=y+3)

## now fit the model...
b <- gam(list(y~s(x1)+s(x2),~s(x1)+s(x2)),family=multinom(K=2))

# Try with weights
reg_weights <- sample(1:100, length(y),replace=T)
b_2 <- gam(list(y~s(x1)+s(x2),~s(x1)+s(x2)),family=multinom(K=2), weights = 
reg_weights)
b; b_2 #identical

1 Ответ

0 голосов
/ 11 мая 2018

Обновление - я думаю, что пакет VGAM будет делать то, что я хочу. Вот рабочий пример использования данных пневмо от VGAM для создания% в каждом классе. 

library(mgcv)
library(VGAM)

# Get pneumo data and log transform exposure time
pneumo <- transform(pneumo, let = log(exposure.time))

# Generate % in each class - note, not suggesting you do 
# this if you have counts, just as an illustrative example if 
# you have % in each class
props <- round(prop.table(as.matrix(pneumo[,c("normal", "mild", 
"severe")]),margin=1)*100,0)

# Fit vgam mod with smoothed relationship with log exposure time
fit <- vgam(props ~ s(let), propodds, data = pneumo)
summary(fit)
plot(fit)
...