在使用mgcv包运行gam模型时,我遇到了一条奇怪的错误消息,我无法理解:
"model.frame.default中的错误(公式=死亡~pm10 +滞后(resid1,1)+:变量长度不同(找到'Lag(resid1,1)')".
模型1中使用的观察数量与偏差残差的长度完全相同,因此我认为此误差与数据大小或长度的差异无关.
我在网上找到了一个相当有关的错误信息在这里,但后没有得到充分的答案,因此它是不利于我的问题.
可重复的示例和数据如下:
library(quantmod)
library(mgcv)
require(dlnm)
df <- chicagoNMMAPS
df1 <- df[,c("date","dow","death","temp","pm10")]
df1$trend<-seq(dim(df1)[1]) ### Create a time trend
model1<-gam(death ~ pm10 + s(trend,k=14*7)+ s(temp,k=5),
data=df1, na.action=na.omit, family=poisson)
resid1 <- residuals(model1,type="deviance")
model1_1 <- update(model1,.~.+ Lag(resid1,1), na.action=na.omit)
model1_2<-gam(death ~ pm10 + s(trend,k=14*7)+ s(temp,k=5) + Lag(resid1,1), data=df1,
na.action=na.omit, family=poisson)
这两个模型都产生了相同的错误消息.
我有几年的时间序列,我需要在一个图表中绘制.最大的系列平均值为340,最小值为245,最大值为900.最小系列的平均值为7,最小值为-28,最大值为31.其余系列的值均为6到700.该系列经历了多年的常规年度和季节性模式,直到突然出现一个月的温度升高,随后死亡人数大大增加.
我无法提供任何实际数据,但我已经模拟了以下数据并尝试了下面的代码,该代码基于此处的示例代码http://www.r-bloggers.com/multiple-y-axis-in-ar-情节/.但情节并没有产生我想要的东西.我有以下问题
非常感谢
temp<- rnorm(365, 5, 10)
mort<- rnorm(365, 300, 45)
poll<- rpois(365, lambda=76)
date<-seq(as.Date('2011-01-01'),as.Date('2011-12-31'),by = 1)
df<-data.frame(date,mort,poll,temp)
windows(600,600)
par(mar=c(5, 12, 4, 4) + 0.1)
with(df, {
plot(date, mort, axes=F, ylim=c(170,max(mort)), xlab="", ylab="",type="l",col="black", main="")
points(date,mort,pch=20,col="black")
axis(2, ylim=c(170,max(mort)),col="black",lwd=2)
mtext(2,text="Mortality",line=2)
})
par(new=T)
plot(date, poll, axes=F, ylim=c(45,max(poll)), xlab="", ylab="",
type="l",col="red",lty=2, main="",lwd=1)
axis(2, ylim=c(45,max(poll)),lwd=1,line=3.5)
points(date, poll,pch=20)
mtext(2,text="PM10",line=5.5)
par(new=T)
plot(date, temp, axes=F, ylim=c(-28,max(temp)), xlab="", ylab="",
type="l",lty=3,col="brown", main="",lwd=1)
axis(2, ylim=c(-28,max(temp)),lwd=1,line=7)
points(date, temp,pch=20)
mtext(2,text="Temperature",line=9)
axis(1,pretty(range(date),10))
mtext("date",side=1,col="black",line=2)
我很难将图例文本对齐到左侧.
library(ggplot2)
library(reshape2)
o3<- rnorm(1827, 50, 10)
NO2<- rnorm(1827, 35, 7)
NOx<- rnorm(1827, 45, 10)
pm25<- rnorm(1827, 15, 4)
date<-seq(as.Date('2000-01-01'),as.Date('2004-12-31'),by = 1)
df<-data.frame(date,o3,NO2,NOx,pm25)
meltdf <- melt(df,id="date")
使用此代码,对齐自动向左
ggplot(meltdf, aes(x = date, y = value, colour =variable)) + geom_smooth() + stat_smooth(method = "gam")
然而,以下是alignemt到中心.
ggplot(meltdf, aes(x = date, y = value, colour =variable)) +
geom_smooth() + stat_smooth(method = "gam") +
scale_color_discrete(name="Pollutant" ,labels = c(expression(O[3]),
expression(NO[2]),
expression(NO[x]),
expression(PM[2.5])))
我怎么能用最后一个脚本实现左对齐？
使用mgcv的惩罚样条,我希望在示例数据中获得10 /年的有效自由度(EDF)(整个期间为60).
library(mgcv)
library(dlnm)
df <- chicagoNMMAPS
df1<-subset(df, as.Date(date) >= '1995-01-01')
mod1 <-gam(resp ~ s(time,bs='cr',k=6*15, fx=F)+ s(temp,k=6, bs='cr') + as.factor(dow)
,family=quasipoisson,na.action=na.omit,data=df1)
在示例数据中,由edf测量的时间的基本维度为56.117,小于每年10.
summary(mod1)
Approximate significance of smooth terms:
edf Ref.df F p-value
s(time) 56.117 67.187 5.369 <2e-16 ***
s(temp) 2.564 3.204 0.998 0.393
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.277 Deviance explained = 28.2%
GCV score = 1.1297 Scale est. = 1.0959 n = 2192
手动我将通过提供如下的平滑参数来更改edf a
mod1$sp
s(time) …我是R的新手,目前正在读一本书"广义的添加模型",一本R by Wood(2006)的介绍,并经历了一些练习,特别是关于空气污染和死亡的部分,这是我感兴趣的领域.使用mgcv包我运行以下模型.
library(gamair)
library(mgcv)
data(chicago)
ap1<-gam(death ~ pm10median + so2median + o3median +s(time,bs="cr",k=200)+ s(tmpd,bs="cr"), data=chicago,family=poisson)
如何提取pm10median的效果估计值和x的95%CI并将输出导出为CSV或任何其他选项？
我估计了比值比,相应的95%CI的六种污染物超过4个滞后期.如何创建一个类似于R中附图的垂直图？下图是在SPSS中创建的.生成该图的示例数据如下:
lag pollut or lcl ucl
0 CO 0.97 0.90 1.06
0 PM10 1.00 0.91 1.09
0 NO 0.97 0.92 1.02
0 NO2 1.01 0.89 1.15
0 SO2 0.97 0.85 1.11
0 Ozone 1.00 0.87 1.15
1 CO 1.03 0.95 1.10
1 PM10 0.93 0.86 1.01
1 NO 1.01 0.97 1.06
1 NO2 1.08 0.97 1.20
1 SO2 0.94 0.84 1.04
1 Ozone 0.94 0.84 1.04
2 CO 1.09 1.02 1.16
2 PM10 1.04 0.96 1.13
2 NO 1.04 1.00 …我有多个重复的ID - 日期集,并希望删除除第一个之外的所有重复集.我怎么能这样做？论坛上没有一个例子为我的案例提供了一个有效的例子.根据提供的输出数据,我希望保留以下内容:
1 767 10-dec-97 1
2 767 10-dec-97 2
3 767 10-dec-97 3
4 767 10-dec-97 4
9 19025 11-dec-97 1
10 19025 11-dec-97 2
11 19025 11-dec-97 3
18 27452 16-apr-95 1
19 27452 16-apr-95 2
20 27452 16-apr-95 3
21 27452 16-apr-95 4
示例数据如下:
structure(list(id = c(767L, 767L, 767L, 767L, 9271L, 9271L, 9271L,
9271L, 19025L, 19025L, 19025L, 162749L, 162749L, 162749L, 183446L,
183446L, 183446L, 27452L, 27452L, 27452L, 27452L, 84002L, 84002L,
84002L, 84002L, 276172L, 276172L, 276172L, 276172L), date …