[R] gam and optim

[Greg Dropkin]

hi

probably a silly mistake, but I expected gam to minimise the penalised
deviance.

thanks

greg

set.seed(1)
library(mgcv)
x<-runif(100)
lp<-exp(-2*x)*sin(8*x)
y<-rpois(100,exp(lp))
plot(x,y)
m1<-gam(y~s(x),poisson)
points(x,exp(lp),pch=16,col="green3")
points(x,fitted(m1),pch=16,cex=0.5,col="blue")
W<-diag(fitted(m1))
X<-predict(m1,type="lpmatrix")
S<-m1$smooth[[1]]$S[[1]]
S<-rbind(0,cbind(0,S))
A<-X%*%solve(t(X)%*%W%*%X+m1$sp*S)%*%t(X)%*%W
sum(diag(A))
sum(m1$edf)
fit<-fitted(m1)
b<-m1$coef
range(exp(X%*%b)-fit)
z<-y/fit-1+X%*%b
range(A%*%z-X%*%b)

dv<-function(t)
{
f<-exp(X%*%t)
-2*sum(y*log(f)-f-ifelse(y==0,0,y*log(y))+y)+t%*%S%*%t
}
dv(b)
m1$dev+b%*%S%*%b

#so far, so good


t1<-optim(rep(0,10),dv)
t1$p
b

#different

dv(t1$p)
dv(b)

#different, and dv(t1$p) is lower!

fit1<-exp(X%*%t1$p)
points(x,fit1,pch=16,cex=0.5,col="red")

# different
# gam found b which does approximate the true curve, but does not minimise
the penalised deviance, by a long shot.