[R] graph: add 2 inches on the left outer region, but keep everything unchanged

2014-04-17 Thread Xing Zhao
Hi R experts,

My original graph was plotted, and for some reason, I need to add
extra '2' inches on the left side.
Meanwhile, I want to keep everything unchanged. Particularly, the
length-width ratio for each panel of the original graph is nice,
therefore I want to keep the original ratio

Adding 2 inches to the pdf(width=) and oma=c(0,2,0,0) does not keep
the original length-width ratio.


Thanks for your help
Xing


#orignal plot

pdf(file="d:/test.pdf",width=7, height=7)
par(mfrow=c(3,4), mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.037)

plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')

dev.off()


#new plot
#want to keep everything unchanged, but 2 inches on the left outer region

pdf(file="d:/test.pdf",width=9, height=7)
par(mfrow=c(3,4), mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.03,oma=c(0,2,0,0))

plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')
plot(c(-2,32),c(-0.1,0.9),
type="line",ylim=c(-0.1,0.9),xlab='',ylab="", xaxt='n', yaxt='n')

dev.off()

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] The explanation of ns() with df =2

2014-04-16 Thread Xing Zhao
Dear John,

Thanks for your patience, got your idea

Best,
Xing

On Wed, Apr 16, 2014 at 4:11 AM, John Fox  wrote:
> Dear Xing,
>
> My point wasn't that you should add 1 df for a constant to each cubic but 
> that you need not subtract 1 for continuity. I'm sorry that I didn't make 
> that clear. When you use the natural spline in a regression, there's a 
> constant in the model as a separate term with 1 df, which in effect adjusts 
> the level of the regression curve (assuming one X for simplicity). If you add 
> a constant to each component cubic, the constant in the model is rendered 
> redundant, accounting for the extra df.
>
> I'm afraid that if this is still not sufficiently clear, I'll have to defer 
> to someone with greater powers of explanation.
>
> Best,
> John
>
> -
> John Fox
> McMaster University
> Hamilton, Ontario, Canada
> http://socserv.socsci.mcmaster.ca/jfox
>
>> On Apr 16, 2014, at 12:10 AM, Xing Zhao  wrote:
>>
>> Dear John
>>
>> Sorry I use 3 degree of freedom for  cubic spline. After using 4, it
>> is still not 2. I may make some naive mistake, but I cannot figure
>> out. Where is the problem?
>>
>> 4 (cubic on the right side of the *interior* knot 8)
>> + 4 (cubic on the left side of the *interior* knot 8)
>> - 1 (two curves must be continuous at the *interior* knot 8)
>> - 1 (two curves must have 1st order derivative continuous at the
>> *interior* knot 8)
>> - 1 (two curves must have 2nd order derivative continuous at the
>> *interior* knot 8)
>> - 1 (right side cubic curve must have 2nd order derivative = 0 at the
>> boundary knot 15 due to the linearity constraint)
>> - 1 (similar for the left)
>> = 3, not 2
>>
>> Thanks
>> Xing
>>
>>> On Tue, Apr 15, 2014 at 10:54 AM, John Fox  wrote:
>>> Dear Xing,
>>>
>>>> -Original Message-
>>>> From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
>>>> project.org] On Behalf Of Xing Zhao
>>>> Sent: Tuesday, April 15, 2014 1:18 PM
>>>> To: John Fox
>>>> Cc: r-help@r-project.org; Michael Friendly
>>>> Subject: Re: [R] The explanation of ns() with df =2
>>>>
>>>> Dear Michael and Fox
>>>>
>>>> Thanks for your elaboration. Combining your explanations would, to my
>>>> understanding, lead to the following  calculation of degree of
>>>> freedoms.
>>>>
>>>> 3 (cubic on the right side of the *interior* knot 8)
>>>> + 3 (cubic on the left side of the *interior* knot 8)
>>>> - 1 (two curves must be continuous at the *interior* knot 8)
>>>
>>> You shouldn't subtract 1 for continuity since you haven't allowed a
>>> different level on each side of the knot (that is your initial counting of 3
>>> parameters for the cubic doesn't include a constant).
>>>
>>> Best,
>>> John
>>>
>>>> - 1 (two curves must have 1st order derivative continuous at the
>>>> *interior* knot 8)
>>>> - 1 (two curves must have 2nd order derivative continuous at the
>>>> *interior* knot 8)
>>>> - 1 (right side cubic curve must have 2nd order derivative = 0 at the
>>>> boundary knot 15 due to the linearity constraint)
>>>> - 1 (similar for the left)
>>>> = 1, not 2
>>>>
>>>> Where is the problem?
>>>>
>>>> Best,
>>>> Xing
>>>>
>>>>> On Tue, Apr 15, 2014 at 6:17 AM, John Fox  wrote:
>>>>> Dear Xing Zhao,
>>>>>
>>>>> To elaborate slightly on Michael's comments, a natural cubic spline
>>>> with 2 df has one *interior* knot and two boundary knots (as is
>>>> apparent in the output you provided). The linearity constraint applies
>>>> beyond the boundary knots.
>>>>>
>>>>> I hope this helps,
>>>>> John
>>>>>
>>>>> 
>>>>> John Fox, Professor
>>>>> McMaster University
>>>>> Hamilton, Ontario, Canada
>>>>> http://socserv.mcmaster.ca/jfox/
>>>>>
>>>>> On Tue, 15 Apr 2014 08:18:40 -0400
>>>>> Michael Friendly  wrote:
>>>>>> No, the curves on each side of the know are cubics, joined
>>>>>> so they are continuous.  Se the discussion in \S 17.2 in
>>&g

Re: [R] The explanation of ns() with df =2

2014-04-15 Thread Xing Zhao
Dear John

Sorry I use 3 degree of freedom for  cubic spline. After using 4, it
is still not 2. I may make some naive mistake, but I cannot figure
out. Where is the problem?

4 (cubic on the right side of the *interior* knot 8)
+ 4 (cubic on the left side of the *interior* knot 8)
- 1 (two curves must be continuous at the *interior* knot 8)
- 1 (two curves must have 1st order derivative continuous at the
*interior* knot 8)
- 1 (two curves must have 2nd order derivative continuous at the
*interior* knot 8)
- 1 (right side cubic curve must have 2nd order derivative = 0 at the
boundary knot 15 due to the linearity constraint)
- 1 (similar for the left)
= 3, not 2

Thanks
Xing

On Tue, Apr 15, 2014 at 10:54 AM, John Fox  wrote:
> Dear Xing,
>
>> -Original Message-
>> From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
>> project.org] On Behalf Of Xing Zhao
>> Sent: Tuesday, April 15, 2014 1:18 PM
>> To: John Fox
>> Cc: r-help@r-project.org; Michael Friendly
>> Subject: Re: [R] The explanation of ns() with df =2
>>
>> Dear Michael and Fox
>>
>> Thanks for your elaboration. Combining your explanations would, to my
>> understanding, lead to the following  calculation of degree of
>> freedoms.
>>
>> 3 (cubic on the right side of the *interior* knot 8)
>> + 3 (cubic on the left side of the *interior* knot 8)
>> - 1 (two curves must be continuous at the *interior* knot 8)
>
> You shouldn't subtract 1 for continuity since you haven't allowed a
> different level on each side of the knot (that is your initial counting of 3
> parameters for the cubic doesn't include a constant).
>
> Best,
>  John
>
>> - 1 (two curves must have 1st order derivative continuous at the
>> *interior* knot 8)
>> - 1 (two curves must have 2nd order derivative continuous at the
>> *interior* knot 8)
>> - 1 (right side cubic curve must have 2nd order derivative = 0 at the
>> boundary knot 15 due to the linearity constraint)
>> - 1 (similar for the left)
>> = 1, not 2
>>
>> Where is the problem?
>>
>> Best,
>> Xing
>>
>> On Tue, Apr 15, 2014 at 6:17 AM, John Fox  wrote:
>> > Dear Xing Zhao,
>> >
>> > To elaborate slightly on Michael's comments, a natural cubic spline
>> with 2 df has one *interior* knot and two boundary knots (as is
>> apparent in the output you provided). The linearity constraint applies
>> beyond the boundary knots.
>> >
>> > I hope this helps,
>> >  John
>> >
>> > 
>> > John Fox, Professor
>> > McMaster University
>> > Hamilton, Ontario, Canada
>> > http://socserv.mcmaster.ca/jfox/
>> >
>> > On Tue, 15 Apr 2014 08:18:40 -0400
>> >  Michael Friendly  wrote:
>> >> No, the curves on each side of the know are cubics, joined
>> >> so they are continuous.  Se the discussion in \S 17.2 in
>> >> Fox's Applied Regression Analysis.
>> >>
>> >> On 4/15/2014 4:14 AM, Xing Zhao wrote:
>> >> > Dear all
>> >> >
>> >> > I understand the definition of Natural Cubic Splines are those
>> with
>> >> > linear constraints on the end points. However, it is hard to think
>> >> > about how this can be implement when df=2. df=2 implies there is
>> just
>> >> > one knot, which, according the the definition, the curves on its
>> left
>> >> > and its right should be both be lines. This means the whole line
>> >> > should be a line. But when making some fits. the result still
>> looks
>> >> > like 2nd order polynomial.
>> >> >
>> >> > How to think about this problem?
>> >> >
>> >> > Thanks
>> >> > Xing
>> >> >
>> >> > ns(1:15,df =2)
>> >> >1   2
>> >> >   [1,] 0.000  0.
>> >> >   [2,] 0.1084782 -0.07183290
>> >> >   [3,] 0.2135085 -0.13845171
>> >> >   [4,] 0.3116429 -0.19464237
>> >> >   [5,] 0.3994334 -0.23519080
>> >> >   [6,] 0.4734322 -0.25488292
>> >> >   [7,] 0.5301914 -0.24850464
>> >> >   [8,] 0.5662628 -0.21084190
>> >> >   [9,] 0.5793481 -0.13841863
>> >> > [10,] 0.5717456 -0.03471090
>> >> > [11,] 0.5469035  0.09506722
>> >> > [12,] 0.5082697  0.24570166
>> >> > [13,] 0.4592920  0.41197833
>

Re: [R] The explanation of ns() with df =2

2014-04-15 Thread Xing Zhao
Dear Michael and Fox

Thanks for your elaboration. Combining your explanations would, to my
understanding, lead to the following  calculation of degree of
freedoms.

3 (cubic on the right side of the *interior* knot 8)
+ 3 (cubic on the left side of the *interior* knot 8)
- 1 (two curves must be continuous at the *interior* knot 8)
- 1 (two curves must have 1st order derivative continuous at the
*interior* knot 8)
- 1 (two curves must have 2nd order derivative continuous at the
*interior* knot 8)
- 1 (right side cubic curve must have 2nd order derivative = 0 at the
boundary knot 15 due to the linearity constraint)
- 1 (similar for the left)
= 1, not 2

Where is the problem?

Best,
Xing

On Tue, Apr 15, 2014 at 6:17 AM, John Fox  wrote:
> Dear Xing Zhao,
>
> To elaborate slightly on Michael's comments, a natural cubic spline with 2 df 
> has one *interior* knot and two boundary knots (as is apparent in the output 
> you provided). The linearity constraint applies beyond the boundary knots.
>
> I hope this helps,
>  John
>
> 
> John Fox, Professor
> McMaster University
> Hamilton, Ontario, Canada
> http://socserv.mcmaster.ca/jfox/
>
> On Tue, 15 Apr 2014 08:18:40 -0400
>  Michael Friendly  wrote:
>> No, the curves on each side of the know are cubics, joined
>> so they are continuous.  Se the discussion in \S 17.2 in
>> Fox's Applied Regression Analysis.
>>
>> On 4/15/2014 4:14 AM, Xing Zhao wrote:
>> > Dear all
>> >
>> > I understand the definition of Natural Cubic Splines are those with
>> > linear constraints on the end points. However, it is hard to think
>> > about how this can be implement when df=2. df=2 implies there is just
>> > one knot, which, according the the definition, the curves on its left
>> > and its right should be both be lines. This means the whole line
>> > should be a line. But when making some fits. the result still looks
>> > like 2nd order polynomial.
>> >
>> > How to think about this problem?
>> >
>> > Thanks
>> > Xing
>> >
>> > ns(1:15,df =2)
>> >1   2
>> >   [1,] 0.000  0.
>> >   [2,] 0.1084782 -0.07183290
>> >   [3,] 0.2135085 -0.13845171
>> >   [4,] 0.3116429 -0.19464237
>> >   [5,] 0.3994334 -0.23519080
>> >   [6,] 0.4734322 -0.25488292
>> >   [7,] 0.5301914 -0.24850464
>> >   [8,] 0.5662628 -0.21084190
>> >   [9,] 0.5793481 -0.13841863
>> > [10,] 0.5717456 -0.03471090
>> > [11,] 0.5469035  0.09506722
>> > [12,] 0.5082697  0.24570166
>> > [13,] 0.4592920  0.41197833
>> > [14,] 0.4034184  0.58868315
>> > [15,] 0.3440969  0.77060206
>> > attr(,"degree")
>> > [1] 3
>> > attr(,"knots")
>> > 50%
>> >8
>> > attr(,"Boundary.knots")
>> > [1]  1 15
>> > attr(,"intercept")
>> > [1] FALSE
>> > attr(,"class")
>> > [1] "ns" "basis"  "matrix"
>> >
>>
>>
>> --
>> Michael Friendly Email: friendly AT yorku DOT ca
>> Professor, Psychology Dept. & Chair, Quantitative Methods
>> York University  Voice: 416 736-2100 x66249 Fax: 416 736-5814
>> 4700 Keele StreetWeb:   http://www.datavis.ca
>> Toronto, ONT  M3J 1P3 CANADA
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>
>
>
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] The explanation of ns() with df =2

2014-04-15 Thread Xing Zhao
Dear all

I understand the definition of Natural Cubic Splines are those with
linear constraints on the end points. However, it is hard to think
about how this can be implement when df=2. df=2 implies there is just
one knot, which, according the the definition, the curves on its left
and its right should be both be lines. This means the whole line
should be a line. But when making some fits. the result still looks
like 2nd order polynomial.

How to think about this problem?

Thanks
Xing

ns(1:15,df =2)
  1   2
 [1,] 0.000  0.
 [2,] 0.1084782 -0.07183290
 [3,] 0.2135085 -0.13845171
 [4,] 0.3116429 -0.19464237
 [5,] 0.3994334 -0.23519080
 [6,] 0.4734322 -0.25488292
 [7,] 0.5301914 -0.24850464
 [8,] 0.5662628 -0.21084190
 [9,] 0.5793481 -0.13841863
[10,] 0.5717456 -0.03471090
[11,] 0.5469035  0.09506722
[12,] 0.5082697  0.24570166
[13,] 0.4592920  0.41197833
[14,] 0.4034184  0.58868315
[15,] 0.3440969  0.77060206
attr(,"degree")
[1] 3
attr(,"knots")
50%
  8
attr(,"Boundary.knots")
[1]  1 15
attr(,"intercept")
[1] FALSE
attr(,"class")
[1] "ns" "basis"  "matrix"

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] mgcv, should include a intercept for the 'by' varying coefficient model, which is unconstrained

2014-03-17 Thread Xing Zhao
Dear Dr. Wood and other mgcv experts


In ?gam.models, it says that the numeric "by" variable is genrally not
subjected to an identifiability constraint, and I used the example in
?gam.models, finding some differences (code below).

I think the the problem might become serious when several varying
coefficient terms are specified in one model, such as gam(y ~
s(x0,by=x1) + s(x0,by=x2) + s(x0,by=x3),data=dat). In this case, those
three terms are all not constraint, as they generally will not meet
the three conditions for constraint.

I can still implement it like gam(y ~ s(x0,by=x1) + s(x0,by=x2) +
s(x0,by=x3),data=dat), but is it safe? Is there a best way to
implement the model?

Thank you for your help
Best,
Xing


require(mgcv)
set.seed(10)
## simulate date from y = f(x2)*x1 + error
dat <- gamSim(3,n=400)

b<-gam(y ~ s(x2,by=x1),data=dat)
b1<-gam(y ~ s(x2,by=x1)-1,data=dat)

> range(fitted(b)-fitted(b1))
[1] -0.13027648  0.08117196
> summary(dat$f-fitted(b))
   Min. 1st Qu.  MedianMean 3rd Qu.Max.
-0.5265  0.2628  1.2290  1.7710  2.6280  8.8580
> summary(dat$f-fitted(b1))
   Min. 1st Qu.  MedianMean 3rd Qu.Max.
-0.4618  0.2785  1.2250  1.7390  2.5370  8.7310
> summary(dat$y-fitted(b))
Min.  1st Qu.   Median Mean  3rd Qu. Max.
-6.23500 -1.32700 -0.06752  0.0  1.54900  7.01800
> summary(dat$y-fitted(b1))
Min.  1st Qu.   Median Mean  3rd Qu. Max.
-6.26700 -1.40300 -0.09908 -0.03199  1.51900  6.96700

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] How to know the what functions have used datasets in loaded packages.

2014-03-10 Thread Xing Zhao
Hi all,

Once you start your R program, there are example data sets available
within R along with loaded packages. Command data() will list all the
datasets in loaded packages.

Is there a method to know what R functions have used one of the
datasets to present the function's utility, maybe at least for the
basic R distribution.


Thanks,
Xing

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Safe prediction does not work for bivariate polynomial terms?

2014-01-23 Thread Xing Zhao
Hi everyone,

R documents says the safe prediction is implemented, when basis
functions are used, such as poly(), bs(), ns()

This works for univariate basis, but fails in my bivariate polynomial setting.
Can anyone explain the reason?


The following is a small example.

set.seed(731)
x<-runif(300)
y<-runif(300)

f <- function(x, y) {  0.5+2*x+13*x^2-14*y^2+12*x*y+y }

z <- f(x,y)+rnorm(length(x))*0.2

# default orthogonal polynomials basis
mod <- lm (z ~ poly(x,y,degree = 2))

# raw polynomials basis
mod1 <- lm (z ~ poly(x,y,degree = 2, raw = T))

# data points to evaluate, just the first five data
new <- data.frame(x=x[1:5],y= y[1:5])

z[1:5]
[1]  9.796620 10.397851  1.280832  4.028284  4.811709

# two predicted values differ dramatically, and the orthogonal
polynomials basis fails
predict(mod, new)
1 2 3 4 5
121.46776  40.85002  18.67273  31.82417  20.81673
predict(mod1, new)
1 2 3 4 5
 9.981091 10.418628  1.223148  4.031664  4.837099

# However, the fitted.values are the same
mod$fitted.values[1:5]
1 2 3 4 5
 9.981091 10.418628  1.223148  4.031664  4.837099
mod1$fitted.values[1:5]
1 2 3 4 5
 9.981091 10.418628  1.223148  4.031664  4.837099


Thanks in advance
Xing

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


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

2014-01-18 Thread Xing Zhao
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)  group2group1:xgroup2: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.00
2 -18.098424  0.00
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"))
 groupgroup: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

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


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

2014-01-18 Thread Xing Zhao
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)  group2group1:xgroup2: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.00
2 -18.098424  0.00
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"))
 groupgroup: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

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


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

2014-01-18 Thread Xing Zhao
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)  group2group1:xgroup2: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.00
2 -18.098424  0.00
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"))
 groupgroup: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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.