Exercise 7 Growth Models - Qingzhou Shi
Qingzhou Shi
10/11/2019
Importing Data
library(haven)
tolerance_origin <- read_sav("~/Box Sync/Courses/2019 Fall/BER642/Exercise/7/tolerance.sav")
View(tolerance_origin)
summary(tolerance_origin)## id tol11 tol12 tol13
## Min. : 9.0 Min. :1.000 Min. :1.000 Min. :1.120
## 1st Qu.: 410.0 1st Qu.:1.120 1st Qu.:1.195 1st Qu.:1.310
## Median : 673.5 Median :1.220 Median :1.340 Median :1.725
## Mean : 762.8 Mean :1.364 Mean :1.441 Mean :1.676
## 3rd Qu.:1009.8 3rd Qu.:1.450 3rd Qu.:1.700 3rd Qu.:1.990
## Max. :1653.0 Max. :2.230 Max. :1.990 Max. :2.230
## tol14 tol15 male exposure
## Min. :1.000 Min. :1.120 Min. :0.0000 Min. :0.8100
## 1st Qu.:1.450 1st Qu.:1.310 1st Qu.:0.0000 1st Qu.:0.9225
## Median :1.730 Median :1.945 Median :0.0000 Median :1.1450
## Mean :1.754 Mean :1.861 Mean :0.4375 Mean :1.1912
## 3rd Qu.:1.990 3rd Qu.:2.120 3rd Qu.:1.0000 3rd Qu.:1.3950
## Max. :3.460 Max. :3.320 Max. :1.0000 Max. :1.9900structure(tolerance_origin)## # A tibble: 16 x 8
## id tol11 tol12 tol13 tol14 tol15 male exposure
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 9 2.23 1.79 1.9 2.12 2.66 0 1.54
## 2 45 1.12 1.45 1.45 1.45 1.99 1 1.16
## 3 268 1.45 1.34 1.99 1.79 1.34 1 0.9
## 4 314 1.22 1.22 1.55 1.12 1.12 0 0.81
## 5 442 1.45 1.99 1.45 1.67 1.9 0 1.13
## 6 514 1.34 1.67 2.23 2.12 2.44 1 0.9
## 7 569 1.79 1.9 1.9 1.99 1.99 0 1.99
## 8 624 1.12 1.12 1.22 1.12 1.22 1 0.98
## 9 723 1.22 1.34 1.12 1 1.12 0 0.81
## 10 918 1 1 1.22 1.99 1.22 0 1.21
## 11 949 1.99 1.55 1.12 1.45 1.55 1 0.93
## 12 978 1.22 1.34 2.12 3.46 3.32 1 1.59
## 13 1105 1.34 1.9 1.99 1.9 2.12 1 1.38
## 14 1542 1.22 1.22 1.99 1.79 2.12 0 1.44
## 15 1552 1 1.12 2.23 1.55 1.55 0 1.04
## 16 1653 1.11 1.11 1.34 1.55 2.12 0 1.25Restructuring Data
Restructuring Data from wide format to long format using ‘tolerance.sav’ data file. Save the restructured file using new file name: tolerance_long.sav.
library(reshape)
longFM <- melt(as.data.frame(tolerance_origin), id=c("id","male", "exposure"),Measured=c("tol11", "tol12", "tol13", "tol14", "tol15"))
tolerance_long <- rename(longFM,c("value"="Tolerance","variable"="Time"))
levels(tolerance_long$Time) <- c(1, 2, 3, 4, 5)
structure(tolerance_long)## id male exposure Time Tolerance
## 1 9 0 1.54 1 2.23
## 2 45 1 1.16 1 1.12
## 3 268 1 0.90 1 1.45
## 4 314 0 0.81 1 1.22
## 5 442 0 1.13 1 1.45
## 6 514 1 0.90 1 1.34
## 7 569 0 1.99 1 1.79
## 8 624 1 0.98 1 1.12
## 9 723 0 0.81 1 1.22
## 10 918 0 1.21 1 1.00
## 11 949 1 0.93 1 1.99
## 12 978 1 1.59 1 1.22
## 13 1105 1 1.38 1 1.34
## 14 1542 0 1.44 1 1.22
## 15 1552 0 1.04 1 1.00
## 16 1653 0 1.25 1 1.11
## 17 9 0 1.54 2 1.79
## 18 45 1 1.16 2 1.45
## 19 268 1 0.90 2 1.34
## 20 314 0 0.81 2 1.22
## 21 442 0 1.13 2 1.99
## 22 514 1 0.90 2 1.67
## 23 569 0 1.99 2 1.90
## 24 624 1 0.98 2 1.12
## 25 723 0 0.81 2 1.34
## 26 918 0 1.21 2 1.00
## 27 949 1 0.93 2 1.55
## 28 978 1 1.59 2 1.34
## 29 1105 1 1.38 2 1.90
## 30 1542 0 1.44 2 1.22
## 31 1552 0 1.04 2 1.12
## 32 1653 0 1.25 2 1.11
## 33 9 0 1.54 3 1.90
## 34 45 1 1.16 3 1.45
## 35 268 1 0.90 3 1.99
## 36 314 0 0.81 3 1.55
## 37 442 0 1.13 3 1.45
## 38 514 1 0.90 3 2.23
## 39 569 0 1.99 3 1.90
## 40 624 1 0.98 3 1.22
## 41 723 0 0.81 3 1.12
## 42 918 0 1.21 3 1.22
## 43 949 1 0.93 3 1.12
## 44 978 1 1.59 3 2.12
## 45 1105 1 1.38 3 1.99
## 46 1542 0 1.44 3 1.99
## 47 1552 0 1.04 3 2.23
## 48 1653 0 1.25 3 1.34
## 49 9 0 1.54 4 2.12
## 50 45 1 1.16 4 1.45
## 51 268 1 0.90 4 1.79
## 52 314 0 0.81 4 1.12
## 53 442 0 1.13 4 1.67
## 54 514 1 0.90 4 2.12
## 55 569 0 1.99 4 1.99
## 56 624 1 0.98 4 1.12
## 57 723 0 0.81 4 1.00
## 58 918 0 1.21 4 1.99
## 59 949 1 0.93 4 1.45
## 60 978 1 1.59 4 3.46
## 61 1105 1 1.38 4 1.90
## 62 1542 0 1.44 4 1.79
## 63 1552 0 1.04 4 1.55
## 64 1653 0 1.25 4 1.55
## 65 9 0 1.54 5 2.66
## 66 45 1 1.16 5 1.99
## 67 268 1 0.90 5 1.34
## 68 314 0 0.81 5 1.12
## 69 442 0 1.13 5 1.90
## 70 514 1 0.90 5 2.44
## 71 569 0 1.99 5 1.99
## 72 624 1 0.98 5 1.22
## 73 723 0 0.81 5 1.12
## 74 918 0 1.21 5 1.22
## 75 949 1 0.93 5 1.55
## 76 978 1 1.59 5 3.32
## 77 1105 1 1.38 5 2.12
## 78 1542 0 1.44 5 2.12
## 79 1552 0 1.04 5 1.55
## 80 1653 0 1.25 5 2.12Combined From
Model 1 Level-1: \(Tolerance_{ij}=β_{0j}+β_{1j} Time_{ij}+r_{ij}\) Level-2: \(β_{0j}=γ_{00}+u_{0j}\) \(β_{1j}=γ_{10}+u_{1j}\)
Combined form: \(Tolerance_{ij}=γ_{00}+u_{0j}+γ_{10} Time_{ij}+u_{1j} Time_{ij}+r_{ij}\) \(=(γ_{00}+γ_{10} Time_{ij})+(u_{0j}+u_{1j} Time_{ij}+r_{ij})\)
Model 2 Level-1: \(Tolerance_{ij}=β_{0j}+β_{1j} Time_{ij}+β_{2j} Time_{ij}*Time_{ij}+r_{ij}\) Level-2: \(β_{0j}=γ_{00}+u_{0j}\) \(β_{1j}=γ_{10}+u_{1j}\) \(β_{2j}=γ_{20}+u_{2j}\)
Combined form: \(Tolerance_{ij}=γ_{00}+u_{0j}+γ_{10} Time_{ij}+u_{1j} Time_{ij}+γ_{20} Time_{ij}*Time_{ij}+u_{2j} Time_{ij}*Time_{ij}+r_{ij}\) \(=(γ_{00}+γ_{10} Time_{ij}+γ_{20} Time_{ij}*Time_{ij})+(u_{0j}+u_{1j} Time_{ij}+u_{2j} Time_{ij}*Time_{ij}+r_{ij})\)
Growth Modeling Analyses
Please perform two growth modeling analyses using linear and quadratic trends. Please use the below hierarchical forms. Use maximum likelihood (ML) estimation and AR(1): Heterogeneous covariance type.
Scater Plots
tolerance_long$Time <- as.numeric(tolerance_long$Time)
tolerance_long$Tolerance <- as.numeric(tolerance_long$Tolerance)
tolerance_long$id <- as.factor(tolerance_long$id)
library(ggplot2)
ggplot(tolerance_long,aes(x=Time, y=Tolerance))+geom_point()+geom_smooth()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Model 1
library(nlme)
Model1 <-lme(Tolerance ~ Time, data = tolerance_long,
random= ~ Time | id, method = "ML", na.action = na.exclude, correlation = corAR1(0, form = ~Time|id), control =
list(opt="optim"))summary(Model1)## Linear mixed-effects model fit by maximum likelihood
## Data: tolerance_long
## AIC BIC logLik
## 76.04737 92.72155 -31.02368
##
## Random effects:
## Formula: ~Time | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.08432039 (Intr)
## Time 0.11675100 -0.865
## Residual 0.31947892
##
## Correlation Structure: AR(1)
## Formula: ~Time | id
## Parameter estimate(s):
## Phi
## 0.3160568
## Fixed effects: Tolerance ~ Time
## Value Std.Error DF t-value p-value
## (Intercept) 1.2334012 0.09741695 63 12.661053 0.0000
## Time 0.1283004 0.04051488 63 3.166748 0.0024
## Correlation:
## (Intr)
## Time -0.722
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.6140212 -0.5997018 -0.2236364 0.4787729 2.5639608
##
## Number of Observations: 80
## Number of Groups: 16intervals(Model1)## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 1.04117775 1.2334012 1.4256246
## Time 0.04835634 0.1283004 0.2082445
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: id
## lower est. upper
## sd((Intercept)) 0.0002863866 0.08432039 24.8263339
## sd(Time) 0.0599460550 0.11675100 0.2273844
## cor((Intercept),Time) -1.0000000000 -0.86516271 1.0000000
##
## Correlation structure:
## lower est. upper
## Phi -0.05446026 0.3160568 0.6100712
## attr(,"label")
## [1] "Correlation structure:"
##
## Within-group standard error:
## lower est. upper
## 0.2439544 0.3194789 0.4183847VarCorr(Model1)## id = pdLogChol(Time)
## Variance StdDev Corr
## (Intercept) 0.007109928 0.08432039 (Intr)
## Time 0.013630795 0.11675100 -0.865
## Residual 0.102066781 0.31947892intercept <-gls(Tolerance ~ 1, data = tolerance_long, method = "ML", na.action = na.exclude)
randomIntercept <-lme(Tolerance ~ 1, data = tolerance_long,
random = ~1|id, method = "ML", na.action = na.exclude, control =
list(opt="optim"))
timeRI<-update(randomIntercept, .~. + Time)
timeRS<-update(timeRI, random = ~Time|id)
ARModel<-update(timeRS, correlation = corAR1(0, form = ~Time|id))
summary(ARModel)## Linear mixed-effects model fit by maximum likelihood
## Data: tolerance_long
## AIC BIC logLik
## 76.04737 92.72155 -31.02368
##
## Random effects:
## Formula: ~Time | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.08432039 (Intr)
## Time 0.11675100 -0.865
## Residual 0.31947892
##
## Correlation Structure: AR(1)
## Formula: ~Time | id
## Parameter estimate(s):
## Phi
## 0.3160568
## Fixed effects: Tolerance ~ Time
## Value Std.Error DF t-value p-value
## (Intercept) 1.2334012 0.09741695 63 12.661053 0.0000
## Time 0.1283004 0.04051488 63 3.166748 0.0024
## Correlation:
## (Intr)
## Time -0.722
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.6140212 -0.5997018 -0.2236364 0.4787729 2.5639608
##
## Number of Observations: 80
## Number of Groups: 16Interpretation of Model 1
The resulting model summary and output shows the 95% con- fidence intervals. The effect of time, \(γ_{10} = .13 (.05, .21), t(64) = 3.17, p < .05\), was significant, indicating that tolerance significantly changed over the time period.
The final combined form of Model 1 is as follows: \(Tolerance_{ij}=γ_{00}+u_{0j}+γ_{10} Time_{ij}+u_{1j} Time_{ij}+r_{ij}\) \(=1.23+.13 Time_{ij}\)
Model 2
Model2<-update(Model1, .~. + I(Time^2))
summary(Model2)## Linear mixed-effects model fit by maximum likelihood
## Data: tolerance_long
## AIC BIC logLik
## 77.98368 97.03989 -30.99184
##
## Random effects:
## Formula: ~Time | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.08604446 (Intr)
## Time 0.11692884 -0.855
## Residual 0.31888712
##
## Correlation Structure: AR(1)
## Formula: ~Time | id
## Parameter estimate(s):
## Phi
## 0.3138076
## Fixed effects: Tolerance ~ Time + I(Time^2)
## Value Std.Error DF t-value p-value
## (Intercept) 1.1993605 0.16790518 62 7.143082 0.0000
## Time 0.1588048 0.12893057 62 1.231708 0.2227
## I(Time^2) -0.0050811 0.02038573 62 -0.249246 0.8040
## Correlation:
## (Intr) Time
## Time -0.904
## I(Time^2) 0.812 -0.949
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.6380605 -0.6195373 -0.2166879 0.4655835 2.5610449
##
## Number of Observations: 80
## Number of Groups: 16intervals(Model2)## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 0.87007613 1.199360470 1.52864481
## Time -0.09404516 0.158804808 0.41165477
## I(Time^2) -0.04506018 -0.005081053 0.03489807
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: id
## lower est. upper
## sd((Intercept)) 0.0002249692 0.08604446 32.9096137
## sd(Time) 0.0596816434 0.11692884 0.2290881
## cor((Intercept),Time) -1.0000000000 -0.85491002 1.0000000
##
## Correlation structure:
## lower est. upper
## Phi -0.06723476 0.3138076 0.6149624
## attr(,"label")
## [1] "Correlation structure:"
##
## Within-group standard error:
## lower est. upper
## 0.2413964 0.3188871 0.4212531VarCorr(Model2)## id = pdLogChol(Time)
## Variance StdDev Corr
## (Intercept) 0.007403648 0.08604446 (Intr)
## Time 0.013672353 0.11692884 -0.855
## Residual 0.101688993 0.31888712Interpretation of Model 2
The resulting model summary and output shows the 95% con- fidence intervals. The effect of time, \(γ_{10} = .16 (-.09, .41), t(64) = 1.23, p = .22\), the effect of \(time^2, γ_{20} = .02 (-.05, .03), t(64)=-.25, p=.80\), was both non-significant, indicating that this is not a very good model.
Compare Models
anova(Model1,Model2)## Model df AIC BIC logLik Test L.Ratio p-value
## Model1 1 7 76.04737 92.72155 -31.02368
## Model2 2 8 77.98368 97.03989 -30.99184 1 vs 2 0.06369064 0.8008Anova test between the two models produced a p value of .80 and a very low likelihood ratio of .06, indicating Model 2 is not significantly different from Model 1. Based on all above results, Model 1 is a better fit of the dataset.