Displaying 2 results from an estimated 2 matches for "ldsq".
Did you mean:
lds
2012 Sep 25
1
REML - quasipoisson
...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-fitted(m3))^2/fitted(m3))
phip<-P/(500-sum(m3$edf))
F1p<-Dpq/(2*phip)
F3p<-F3q
#thi...
2012 Oct 01
0
[Fwd: REML - quasipoisson]
...%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-fitted(m3))^2/fitted(m3))
>...