[R-sig-ME] another gls question

[Thierry Onkelinx]

Dear Mark,

SSTOT = SSREG + SSTOT is relation that holds with a lineair model with OLS.
I'm not sure if it still holds with gls.

Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

2015-05-07 22:17 GMT+02:00 Mark Leeds <markleeds2 op gmail.com>:

> Hi All. I'm sorry to bother this list again but I was able to extract the
> orthogonal residuals
> based on what Thierry suggested.  ( Thanks Thierry ).
>
> If I understand correctly, the normalized residuals are the ones that come
> out of
> the transformed regression with independent error term which means that the
> sums of square decomposition should hold when I use them.
>
> Then. since the sums of squares decomposition should hold,  I would expect
> to get the same values for r2= 1-(sse/sstot) and R2 = corr(predicted,
> actual)^2.
>
> Unfortunately, when I use these residuals, the SSTOT = SSREG + SSTOT
> relation still
> isn't holding. ( so of course the r^2 = R^2 relation isn't holding ).
>
>  I've dputted the objects needed to run the code and the code below. Could
> someone take a look and see if I'm doing anything dumb because,
> theoretically, the decomposition should hold. It also could be that I'm
> still not understanding the term "normalized" and what it implies about
> where those residuals are coming from.
>
> My other thought is that maybe ybar need to be modified by subtracting  the
> theta*e_t-1 term
> from it  since the residual is e_t now  rather than e_t + theta* e_t-1 ?
>
> Thanks a lot.
>
> #=================================================================
>
>  dput(x)
> structure(c(3200L, 3084L, 2864L, 2803L, 2755L, 2696L, 2646L,
> 2701L, 2654L, 2766L, 2832L, 2964L, 3041L, 3010L, 3018L, 3374L,
> 3545L, 3441L, 3456L, 3455L, 3503L, 3641L, 3721L, 3828L, 3831L,
> 3858L, 3925L, 3880L, 3935L, 3895L, 3840L, 3756L, 3669L, 3502L,
> 3145L, 2812L, 2586L, 2441L, 234L, 234L, 235L, 237L, 238L, 240L,
> 241L, 242L, 244L, 245L, 246L, 268L, 333L, 335L, 331L, 253L, 243L,
> 241L, 242L, 237L, 242L, 240L, 233L, 232L, 236L, 245L, 256L, 261L,
> 265L, 278L, 291L, 290L, 290L, 307L, 313L, 325L, 339L, 338L), .Dim = c(38L,
> 2L), .Dimnames = list(NULL, c("ewma.ret", "ewma.imbal")))
> > dput(y)
> c(77.1, 92.9, 98.3, 88.1, 79.4, 91, 100.4, 108.9, 123.6, 157.3,
> 154.3, 143.9, 147.5, 97.3, 76.6, 72.8, 68.9, 66.7, 55, 55, 55,
> 54.7, 47.9, 49, 44.4, 40.2, 43.8, 49.5, 50.4, 59.6, 72.4, 70.6,
> 82, 89.5, 101.3, 116.7, 115.2, 122.9)
>
>
> gls.est <- gls(y~x, correlation=corARMA(q=1), method='ML')
>
> ybar <- gls.est$coefficients[1]
> residualsnorm <- residuals(gls.est, type = "normalized")
>
> gls.est.ssenorm <- sum((residualsnorm)^2)
> gls.est.ssreg <-sum((gls.est$fitted - ybar)^2)
> gls.est.sstot <- sum((y - ybar)^2)
>
> gls.est.rsquaredA <- 1 - (gls.est.ssenorm/gls.est.sstot)
> gls.est.rsquaredB <- cor(y, gls.est$fitted)^2
>
> cat("SSE + SSREG: ", gls.est.ssreg + gls.est.ssenorm, "\n")
> cat("SSTOT: ", gls.est.sstot, "\n")
>
> cat("GLS RSQUARED A : ", gls.est.rsquaredA, "\n")
> cat("GLS RSQUARED B : ", gls.est.rsquaredB, "\n")
>
> gls.sum <- summary(gls.est)
> print(gls.sum)
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]