[R] lm(y ~ group/x ) + predict.lm(...,type="terms")

Xing Zhao zhaoxing at uw.edu
Sat Jan 18 20:54:29 CET 2014


Hi, all

I am trying to figure out the computation result for
predict.lm(...,type="terms")  when the original fitting model has a
nesting term, lm(y ~ group/x ).

A example,

> set.seed(731)
> group <- factor(rep(1:2, 200))
> x <- rnorm(400)
>
> fun1 <- function(x) -3*x+8
> fun2 <- function(x) 15*x-18
>
> y <- (group==1)*fun1(x)+(group==2)*fun2(x) + rnorm(400)
>
> mod1 <- lm(y ~ group/(x-1) ) # without intercetp
> mod2 <- lm(y ~ group/x ) # with intercetp
>
>
> #data to predict
> new <- data.frame(x=rep(0:2,each=2),
+                   group=factor(rep(1:2,3)))
> new
  x group
1 0     1
2 0     2
3 1     1
4 1     2
5 2     1
6 2     2
> coef(mod1) # checking coefficients, both make sense to me.
    group1     group2   group1:x   group2:x
  7.864981 -18.098424  -2.963931  15.051449
> coef(mod2)
(Intercept)      group2    group1:x    group2:x
   7.864981  -25.963405   -2.963931   15.051449
>
> predict(mod1, new,type = c("response")) # two "response" type predictions are the same, make sense to me.
         1          2          3          4          5          6
  7.864981 -18.098424   4.901050  -3.046975   1.937120  12.004474
> predict(mod2, new,type = c("response"))
         1          2          3          4          5          6
  7.864981 -18.098424   4.901050  -3.046975   1.937120  12.004474
>
> predict(mod1, new,type = c("terms")) # make sense to me
       group   group:x
1   7.864981  0.000000
2 -18.098424  0.000000
3   7.864981 -2.963931
4 -18.098424 15.051449
5   7.864981 -5.927861
6 -18.098424 30.102898
attr(,"constant")
[1] 0

# I want to know the computation process for group:x below??? this is
what I am interested in
> predict(mod2, new,type = c("terms"))
     group    group:x
1  12.9817  0.5209069
2 -12.9817  0.5209069
3  12.9817 -2.4430237
4 -12.9817 15.5723560
5  12.9817 -5.4069544
6 -12.9817 30.6238052
attr(,"constant")
[1] -5.637629


Thanks in advance
Xing



More information about the R-help mailing list