# load & see
data <- WWGbook::autism

head(data)
##   age vsae sicdegp childid
## 1   2    6       3       1
## 2   3    7       3       1
## 3   5   18       3       1
## 4   9   25       3       1
## 5  13   27       3       1
## 6   2   17       3       3
str(data)
## 'data.frame':    612 obs. of  4 variables:
##  $ age    : int  2 3 5 9 13 2 3 5 9 13 ...
##  $ vsae   : int  6 7 18 25 27 17 18 12 18 24 ...
##  $ sicdegp: int  3 3 3 3 3 3 3 3 3 3 ...
##  $ childid: int  1 1 1 1 1 3 3 3 3 3 ...
#Convert 123 to LMH.
SICLMH <- with(data, cut(sicdegp, ordered=T, breaks=c(0, 1, 2, 3), labels=c("L", "M", "H")))
#cbind
dta <- cbind(data, SICLMH)
head(dta)
##   age vsae sicdegp childid SICLMH
## 1   2    6       3       1      H
## 2   3    7       3       1      H
## 3   5   18       3       1      H
## 4   9   25       3       1      H
## 5  13   27       3       1      H
## 6   2   17       3       3      H
#Def. centered age
centerage <- round((dta$age-mean(dta$age)),2)
ndta <- cbind(dta, centerage)
head(ndta)
##   age vsae sicdegp childid SICLMH centerage
## 1   2    6       3       1      H     -3.77
## 2   3    7       3       1      H     -2.77
## 3   5   18       3       1      H     -0.77
## 4   9   25       3       1      H      3.23
## 5  13   27       3       1      H      7.23
## 6   2   17       3       3      H     -3.77
#draw plot

library(ggplot2)

gg1 <- ggplot(data=ndta, aes(x=centerage, y=vsae)) + labs(x='Age (in years,centered)', y='VASE score') + 
geom_point(alpha = 0.5) + stat_smooth(data=ndta, formula=y ~ x, method='lm', se=T) + 
  facet_wrap(. ~ SICLMH, ncol=3) +
  geom_line(aes(group = childid), alpha = 0.3) + 
  scale_x_continuous(limits = c(-4, 7.5),
                     breaks = seq(-2.5,5,2.5))

gg1
## Warning: Removed 2 rows containing non-finite values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).

## find the data has missing value,then we must delete missing values

nNAdta <- na.omit(ndta)

gg2 <- ggplot(data=nNAdta, aes(x=centerage, y=vsae)) + labs(x='Age (in years,centered)', y='VASE score') + 
geom_point(alpha = 0.5) + stat_smooth(data=nNAdta, formula=y ~ x, method='lm', se=T) + 
  facet_wrap(. ~ SICLMH, ncol=3) +
  geom_line(aes(group = childid), alpha = 0.3) + 
  scale_x_continuous(limits = c(-4, 7.5),
                     breaks = seq(-2.5,5,2.5))

gg2

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
newdata<-nNAdta %>% mutate(age2 = age - 2) %>% 
  # group data
  group_by(SICLMH, age2) %>%
  # compute mean and standard error for vase
  summarize(vsaemean = mean(vsae), 
            vsaese = sd(vsae) / sqrt(n()))

head(newdata)  
## # A tibble: 6 x 4
## # Groups:   SICLMH [2]
##   SICLMH  age2 vsaemean vsaese
##   <ord>  <dbl>    <dbl>  <dbl>
## 1 L          0     7     0.387
## 2 L          1    12.0   0.914
## 3 L          3    15.0   1.47 
## 4 L          7    25.6   4.74 
## 5 L         11    37.1   6.72 
## 6 M          0     8.67  0.435
##draw ggplot
gg3 <- ggplot(data=newdata) +
  aes(age2, vsaemean, group=SICLMH, shape=SICLMH) +
  geom_errorbar(aes(ymin=vsaemean - vsaese,
                    ymax=vsaemean + vsaese),
                width=.2, size=.3, 
                position=position_dodge(.5)) +
  geom_line(position=position_dodge(.5), 
            aes(linetype=SICLMH))+
  geom_point(position=position_dodge(.5), 
             size=rel(2))+
  scale_shape_manual(values = c(1, 2, 19))+
  labs(x="Age (in year -2)", y="VSAE score")+
  theme(legend.position= c(0.07,0.85),)
 
  
  
  
gg3