作者andrew43 (討厭有好心推文後刪文者)
看板R_Language
標題Re: [問題] plot 3D圖如何標示最大值
時間Thu Sep 3 01:11:21 2015
不知道是不是答非所問,你跑看看。
# 建資料
dt.obs <- data.frame(
x = c(3,5,2,1,2,3,8),
y = c(1,6,4,9,4,2,3),
z = c(3,26,9,9,9,5,25)
)
# 迴歸式
mod <- glm(z ~ x + y, data = dt.obs, family=quasipoisson)
summary(mod)
# 放網格和求預測值
library(plot3D)
grid.lines <- 20
M <- mesh(
seq(1, 9, length.out = grid.lines),
seq(1, 9, length.out = grid.lines)
)
dt.pred <- data.frame(
x = as.vector(M$x),
y = as.vector(M$y)
)
dt.pred$z <- predict(mod, dt.pred, type = "response")
# dt.pred 就是網格各點的資料
# 要找特別的 (x,y,z) 可以從這裡找找看
# 一號圖
scatter3D(dt.obs$x, dt.obs$y, dt.obs$z, surf = list(
x = M$x, y = M$y,
z = matrix(dt.pred$z, nrow = grid.lines),
fit = predict(mod, type="response")))
dev.new()
# 二號圖
x <- M$x
y <- M$y
z <- matrix(dt.pred$z, nrow = grid.lines)
surf3D(x,y,z)
※ 引述《aee36900 (持久戰!!)》之銘言:
: [問題類型]:
:
: 程式諮詢(我想用R 做某件事情,但是我不知道要怎麼用R 寫出來)
:
: [軟體熟悉度]:
: 新手(沒寫過程式,R 是我的第一次)
: [問題敘述]:
: 練習dlnm模型
: 想要在3D圖形上將RR最大值的數值及對應的天數找出來
: 不知道有無function可以直接使用計算出來
: 圖形如下
: http://imgur.com/6IcvqCW
: [程式範例]:
:
: install.packages("dlnm")
: library(dlnm)
: cb3.pm <- crossbasis(chicagoNMMAPS$pm10, lag=1, argvar=list(fun="lin",cen=0),
: arglag=list(fun="strata"))
: varknots <- equalknots(chicagoNMMAPS$temp,fun="bs",df=5,degree=2)
: lagknots <- logknots(30, 3)
: cb3.temp <- crossbasis(chicagoNMMAPS$temp, lag=30, argvar=list(fun="bs",knots=varknots,cen=21), arglag=list(knots=lagknots))
: model3 <- glm(death ~ cb3.pm + cb3.temp + ns(time, 7*14) + dow,family=quasipoisson(), chicagoNMMAPS)
: pred3.temp <- crosspred(cb3.temp, model3, by=1)
: plot(pred3.temp, xlab="Temperature", zlab="RR", theta=200, phi=40, lphi=30,
: main="3D graph of temperature effect")
:
: [環境敘述]:
:
: R version 3.1.3 (2015-03-09)
: Platform: x86_64-redhat-linux-gnu (64-bit)
: Running under: CentOS release 6.5 (Final)
:
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 125.230.67.53
※ 文章網址: https://www.ptt.cc/bbs/R_Language/M.1441213890.A.29C.html
※ 編輯: andrew43 (125.230.67.53), 09/03/2015 01:13:18
→ celestialgod: 一號圖上方要加 dt.pred$z = predict(mod, dt.pred) 09/03 01:47
→ celestialgod: qq..看錯 我自己沒複製到XD 09/03 01:48