Elegant way to express residual calculation in R?

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Elegant way to express residual calculation in R?

John McHenry
    Hi All,

    I am illustrating a simple, two-way ANOVA using the following data and I'm
    having difficulty in expressing the predicted values succinctly in R.

    X<- data.frame(read.table(textConnection("
        Machine.1    Machine.2    Machine.3
        53           61           51
        47           55           51
        46           52           49
        50           58           54
        49           54           50"
    ), header=TRUE))
    rownames(X)<- paste("Operator.", 1:nrow(X), sep="")
    print(X)

    # I'd like to know if there is a more elegant way to calculate the residuals
    # than the following, which seems to be rather a kludge. If you care to read
    # the code you'll see what I mean.

    machine.adjustment<-  colMeans(X) - mean(mean(X))    # length(machine.adjustment)==3
    operator.adjustment<- rowMeans(X) - mean(mean(X))    # length(operator.adjustment)==5
    X.predicted<- numeric(0)
    for (j in 1:ncol(X))
    {
        new.col<- mean(mean(X)) + operator.adjustment + machine.adjustment[j]
        X.predicted<- cbind(X.predicted, new.col)
    }
    print(X.predicted)
    X.residual<- X - X.predicted
    SS.E<- sum( X.residual^2 )

It seems like there ought to be some way of doing that a little bit cleaner ...

Thanks,

Jack.

               
---------------------------------

Bring photos to life! New PhotoMail  makes sharing a breeze.
        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Elegant way to express residual calculation in R?

Liaw, Andy
Not sure if the following would qualify as elegant...  It's easier to work
with a matrix instead of data frame:

x.mat <- data.matrix(X)
xmean <- mean(x.mat)
operator.eff <- ave(x.mat, row(x.mat)) - xmean
machine.eff <- ave(x.mat, col(x.mat)) - xmean
(predicted.x <- xmean + operator.eff + machine.eff)
resid.x <- x.mat - predicted.x
sum(resid.x^2)

You can also use aov(), but you already know that...

Andy

From: John McHenry

>
>     Hi All,
>
>     I am illustrating a simple, two-way ANOVA using the
> following data and I'm
>     having difficulty in expressing the predicted values
> succinctly in R.
>
>     X<- data.frame(read.table(textConnection("
>         Machine.1    Machine.2    Machine.3
>         53           61           51
>         47           55           51
>         46           52           49
>         50           58           54
>         49           54           50"
>     ), header=TRUE))
>     rownames(X)<- paste("Operator.", 1:nrow(X), sep="")
>     print(X)
>
>     # I'd like to know if there is a more elegant way to
> calculate the residuals
>     # than the following, which seems to be rather a kludge.
> If you care to read
>     # the code you'll see what I mean.
>
>     machine.adjustment<-  colMeans(X) - mean(mean(X))    #
> length(machine.adjustment)==3
>     operator.adjustment<- rowMeans(X) - mean(mean(X))    #
> length(operator.adjustment)==5
>     X.predicted<- numeric(0)
>     for (j in 1:ncol(X))
>     {
>         new.col<- mean(mean(X)) + operator.adjustment +
> machine.adjustment[j]
>         X.predicted<- cbind(X.predicted, new.col)
>     }
>     print(X.predicted)
>     X.residual<- X - X.predicted
>     SS.E<- sum( X.residual^2 )
>
> It seems like there ought to be some way of doing that a
> little bit cleaner ...
>
> Thanks,
>
> Jack.
>
>
> ---------------------------------
>
> Bring photos to life! New PhotoMail  makes sharing a breeze.
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Elegant way to express residual calculation in R?

Gabor Grothendieck
In reply to this post by John McHenry
Try this:

Xm <- as.matrix(X)
X.lm <- lm(c(Xm) ~ factor(col(Xm)) + factor(row(Xm)))
sum(resid(X.lm)^2)

or if the idea was to do it without using lm try replacing
your calculation of X.predicted with this:

X.predicted <-
  outer(operator.adjustment, machine.adjustment, "+") + mean(mean(X))


On 2/27/06, John McHenry <[hidden email]> wrote:

>    Hi All,
>
>    I am illustrating a simple, two-way ANOVA using the following data and I'm
>    having difficulty in expressing the predicted values succinctly in R.
>
>    X<- data.frame(read.table(textConnection("
>        Machine.1    Machine.2    Machine.3
>        53           61           51
>        47           55           51
>        46           52           49
>        50           58           54
>        49           54           50"
>    ), header=TRUE))
>    rownames(X)<- paste("Operator.", 1:nrow(X), sep="")
>    print(X)
>
>    # I'd like to know if there is a more elegant way to calculate the residuals
>    # than the following, which seems to be rather a kludge. If you care to read
>    # the code you'll see what I mean.
>
>    machine.adjustment<-  colMeans(X) - mean(mean(X))    # length(machine.adjustment)==3
>    operator.adjustment<- rowMeans(X) - mean(mean(X))    # length(operator.adjustment)==5
>    X.predicted<- numeric(0)
>    for (j in 1:ncol(X))
>    {
>        new.col<- mean(mean(X)) + operator.adjustment + machine.adjustment[j]
>        X.predicted<- cbind(X.predicted, new.col)
>    }
>    print(X.predicted)
>    X.residual<- X - X.predicted
>    SS.E<- sum( X.residual^2 )
>
> It seems like there ought to be some way of doing that a little bit cleaner ...
>
> Thanks,
>
> Jack.
>
>
> ---------------------------------
>
> Bring photos to life! New PhotoMail  makes sharing a breeze.
>        [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Elegant way to express residual calculation in R?

Gabor Grothendieck
Here is one further idea:

Xm <- as.matrix(X)
Xm. <- scale(Xm, scale = FALSE)
Xm.. <- t(scale(t(tmp), scale = FALSE))
sum(Xm..^2)


On 2/27/06, Gabor Grothendieck <[hidden email]> wrote:

> Try this:
>
> Xm <- as.matrix(X)
> X.lm <- lm(c(Xm) ~ factor(col(Xm)) + factor(row(Xm)))
> sum(resid(X.lm)^2)
>
> or if the idea was to do it without using lm try replacing
> your calculation of X.predicted with this:
>
> X.predicted <-
>  outer(operator.adjustment, machine.adjustment, "+") + mean(mean(X))
>
>
> On 2/27/06, John McHenry <[hidden email]> wrote:
> >    Hi All,
> >
> >    I am illustrating a simple, two-way ANOVA using the following data and I'm
> >    having difficulty in expressing the predicted values succinctly in R.
> >
> >    X<- data.frame(read.table(textConnection("
> >        Machine.1    Machine.2    Machine.3
> >        53           61           51
> >        47           55           51
> >        46           52           49
> >        50           58           54
> >        49           54           50"
> >    ), header=TRUE))
> >    rownames(X)<- paste("Operator.", 1:nrow(X), sep="")
> >    print(X)
> >
> >    # I'd like to know if there is a more elegant way to calculate the residuals
> >    # than the following, which seems to be rather a kludge. If you care to read
> >    # the code you'll see what I mean.
> >
> >    machine.adjustment<-  colMeans(X) - mean(mean(X))    # length(machine.adjustment)==3
> >    operator.adjustment<- rowMeans(X) - mean(mean(X))    # length(operator.adjustment)==5
> >    X.predicted<- numeric(0)
> >    for (j in 1:ncol(X))
> >    {
> >        new.col<- mean(mean(X)) + operator.adjustment + machine.adjustment[j]
> >        X.predicted<- cbind(X.predicted, new.col)
> >    }
> >    print(X.predicted)
> >    X.residual<- X - X.predicted
> >    SS.E<- sum( X.residual^2 )
> >
> > It seems like there ought to be some way of doing that a little bit cleaner ...
> >
> > Thanks,
> >
> > Jack.
> >
> >
> > ---------------------------------
> >
> > Bring photos to life! New PhotoMail  makes sharing a breeze.
> >        [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > [hidden email] mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> >
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html