## Function for sampling a new dataset
sampleData <- function(
                       np = 7,
                       mu = 0,
                       sd.p = 4,
                       alpha.p = 8,
                       beta.p = -0.5,
                       pm = rnorm( np, mu.p, sd.p ),
                       ni = sample( 5:10, np, replace=TRUE ),
                       beta.i = 1,
                       sd.i = 1,
                       sd.r = 0.2
                       )
  {
    ## Generating x's in personwise groups around xm
    xm  <- rep( pm, ni )
    pers <- as.integer( factor( xm ) )
    xpi <- xm + rnorm( xm, 0, sd.i )
    
    ## Personwise intercept
    Ap  <- alpha.p + beta.p * xm
    
    ## "Teoretical personwise means"
    ym  <- Ap + beta.i * xm
    ## The model
    y  <-  mu + Ap + beta.i * xpi +  rnorm( xm, 0, sd.r )
    
    data.frame( x=xpi, y=y, p=pers, ym=ym, xm=xm )
  }

db <- sampleData(
                 np=4,
                 pm=c(5,7,9,12),
                 ni=c(10,9,11,12),
                 beta.i=1,
                 beta.p=-4, 
                 alpha.p = 40,
                 sd.i=1.5,
                 sd.r=0.4 )
attach( db , warn.conflicts = FALSE)

## Overall regression
linreg <- lm(y~x)

## Common slope, but groupwise intercept
Mp <- lm( y ~ factor(p) + x - 1 ) 

## Groupwise means
ymeans <- tapply(y, p, mean)

## Scatterplot
plot(x, y,
     ylim = c(-10,40), 
     xlim = c(min(x), max(x)+1),
     xlab = expression(x[1]))

lines(linreg$fitted ~ x, lty = "32")

## Horisontal lines indicating group means
abline(h = ymeans, lty = "11")

  for( i in 1:max(p) )
  {
    ## Regression line
    x1 <- min(x[p==i])
    x2 <- max(x[p==i])
    y2 <- max(y[p==i])
    
    curve(coef(Mp)[i] + coef(Mp)[max(p)+1]*x, x1, x2, add=TRUE)
    text(x2+.8,y2+.8, bquote(x[2]==.(i)),cex=0.7)
  }
detach( db )