Displaying 2 results from an estimated 2 matches for "ldhsq".
Did you mean:
ldhs
2012 Sep 25
1
REML - quasipoisson
...in 2:10)
{
for (j in 2:10)
{
Sq[i,j]=S1q[i-1,j-1]
Sq[i+9,j+9]=S2q[i-1,j-1]
}
}
betaq<-m3$coef
#penalised deviance
Dpq<-m3$dev+t(betaq)%*%Sq%*%betaq
F1q<-Dpq/(2*phiq)
#2nd term undefined
#3rd term
Xq<-predict(m3,type="lpmatrix")
Wq<-diag(fitted(m3))
Hq<-t(Xq)%*%Wq%*%Xq
ldhsq<-determinant(Hq+Sq,log=TRUE)$modulus[1]
eigSq<-eigen(Sq,only.values=TRUE)$val
ldsq<-sum(log(eigSq[1:16]))
F3q<-(ldhsq-ldsq)/2
#4th term
Mp=3
F4q<-Mp/2*log(2*pi*phiq)
F1q+F3q-F4q
m3$gcv
#quite different
#but if phiq is replaced by the Pearson estimate of the scale
P<-sum((y-fit...
2012 Oct 01
0
[Fwd: REML - quasipoisson]
...-1,j-1]
> }
> }
> betaq<-m3$coef
> #penalised deviance
> Dpq<-m3$dev+t(betaq)%*%Sq%*%betaq
> F1q<-Dpq/(2*phiq)
>
> #2nd term undefined
>
> #3rd term
> Xq<-predict(m3,type="lpmatrix")
> Wq<-diag(fitted(m3))
> Hq<-t(Xq)%*%Wq%*%Xq
> ldhsq<-determinant(Hq+Sq,log=TRUE)$modulus[1]
> eigSq<-eigen(Sq,only.values=TRUE)$val
> ldsq<-sum(log(eigSq[1:16]))
> F3q<-(ldhsq-ldsq)/2
>
> #4th term
> Mp=3
> F4q<-Mp/2*log(2*pi*phiq)
>
> F1q+F3q-F4q
> m3$gcv
>
> #quite different
>
> #but if phiq...