Как включить gam() в xyplot() из пакета Lattice?

Я пытаюсь включить обобщенные аддитивные модели gam() из пакета mgcv либо в функцию xyplot(), либо в функцию coplot() из пакета lattice в R.

Данные можно найти в http://statweb.stanford.edu/~tibs/ElemStatLearn/, выбрав данные по озону.

Вот мой код для сглаживания ядра.

    ozonedata=ozone 
    Temperature=equal.count(ozonedata$temperature,4,1/2)
    Wind=equal.count(ozonedata$wind,4,1/2)

    xyplot(ozone^(1/3) ~ radiation | Temperature * Wind, data = ozonedata, as.table = TRUE, 
           panel = function(x, y, ...) {panel.xyplot(x, y, ...);panel.loess(x, y)}, 
           pch = 20,xlab = "Solar Radiation", ylab = "Ozone (ppb)")

or 

    coplot((ozone^(1/3))~radiation|temperature*wind,data=ozonedata,number=c(4,4),
           panel = function(x, y, ...) panel.smooth(x, y, span = .8, ...),
           xlab="Solar radiation (langleys)", ylab="Ozone (cube root ppb)")

Обобщенные аддитивные модели генерируются следующим образом.

    gam_ozone = gam(ozone^(1/3)~s(radiation)+s(temperature)+s(wind),data=ozonedata,method="GCV.Cp")

Теперь у меня возникли проблемы с объединением фитингов из gam() в решетчатые графики.


r gam mgcv lattice
person Yu Deng    schedule 01.06.2014    source источник

Ответы (1)


arrow_upward
4
arrow_downward

Я думаю, это должно сработать. Обратите внимание, что мы передаем объект gam_ozone в качестве параметра объекту xyplot с именем gam. Это позволит нам получить к нему доступ в функции панели.

xyplot(ozone^(1/3) ~ radiation | Temperature * Wind, data = ozonedata, 
    as.table = TRUE, gam=gam_ozone, 
    panel = function(x, y, gam,...) {

        xx <- dimnames(trellis.last.object())
        ww <- arrayInd(packet.number(), .dim=sapply(xx, length))
        avgtmp <- mean(xx[[1]][[ww[1]]])
        avgwind <- mean(xx[[2]][[ww[2]]])

        gx<-seq(min(x, na.rm=T), max(x, na.rm=T), length.out=50)
        gy<-predict(gam, data.frame(radiation=gx, 
            temperature=avgtmp, wind=avgwind))

        panel.xyplot(x, y, ...);
        panel.xyplot(gx, gy, ..., col="red", type="l");
    }, 
    pch = 20,xlab = "Solar Radiation", ylab = "Ozone (ppb)"
)

Теперь, чтобы использовать gam для прогнозирования, вам нужно найти значение, которое будет использоваться для ветра и температуры для каждой панели. Что я решил сделать, так это просто взять среднее значение для каждого диапазона гальки. Поэтому я использовал одобренную Deepayan недокументированную функцию. чтобы получить диапазоны для каждой из текущих черепиц с помощью вызова dimnames. И затем я нахожу текущую панель с помощью packet.number() Когда у меня есть диапазоны, я использую средства для получения среднего значения.

Я не буду использовать модель gam, которую мы передали, чтобы предсказать кривую для каждой панели. Я рассчитываю диапазон из 50 x значений на основе наблюдаемых значений, а затем предсказываю из gam новую строку.

Наконец, я рисую необработанные данные с помощью panel.xyplot, а затем рисую предсказанную линию gam.

пример вывода решетки с плавной кривой гамма

person MrFlick    schedule 01.06.2014