library(effsize)
library(effectsize)## Registered S3 methods overwritten by 'parameters':
## method from
## as.double.parameters_kurtosis datawizard
## as.double.parameters_skewness datawizard
## as.double.parameters_smoothness datawizard
## as.numeric.parameters_kurtosis datawizard
## as.numeric.parameters_skewness datawizard
## as.numeric.parameters_smoothness datawizard
## print.parameters_distribution datawizard
## print.parameters_kurtosis datawizard
## print.parameters_skewness datawizard
## summary.parameters_kurtosis datawizard
## summary.parameters_skewness datawizardlibrary(lmSupport)
library(QuantPsyc)## Loading required package: boot## Loading required package: MASS##
## Attaching package: 'QuantPsyc'## The following object is masked from 'package:base':
##
## normlibrary(plyr)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:boot':
##
## logitlibrary(multicon)## Loading required package: psych##
## Attaching package: 'psych'## The following object is masked from 'package:car':
##
## logit## The following object is masked from 'package:boot':
##
## logit## The following object is masked from 'package:effectsize':
##
## phi## The following object is masked from 'package:effsize':
##
## cohen.d## Loading required package: abind## Loading required package: foreachlibrary(cocor)
library(ggplot2)##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphalibrary(psych)
library(lattice) ##
## Attaching package: 'lattice'## The following object is masked from 'package:boot':
##
## melanomalibrary(plyr)
require(broom)## Loading required package: broomlibrary(nlme)
library(sjmisc)
library(lme4)## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmListsetwd("/Users/ninaradosic/Desktop/Essays/UC Riverside/Longitudinal Analysis/Homework")
HW2_long <- read.csv("HW2_long.csv", na.strings = ".", header=T)
attach(HW2_long)HW2_long$age_centered <- age_centered + 1HW2_cit <- data.frame(id = HW2_long$ID,
var = HW2_long$cit,
age_centered = HW2_long$age_centered,
d_cit = 1,
d_age = 0,
grp = 'cit')
summary(HW2_cit)## id var age_centered d_cit d_age
## Min. : 1.00 Min. :2.111 Min. : 1.00 Min. :1 Min. :0
## 1st Qu.: 50.75 1st Qu.:2.815 1st Qu.: 7.00 1st Qu.:1 1st Qu.:0
## Median :100.50 Median :3.000 Median :14.00 Median :1 Median :0
## Mean :100.50 Mean :3.001 Mean :14.22 Mean :1 Mean :0
## 3rd Qu.:150.25 3rd Qu.:3.185 3rd Qu.:21.00 3rd Qu.:1 3rd Qu.:0
## Max. :200.00 Max. :3.815 Max. :28.00 Max. :1 Max. :0
## grp
## Length:800
## Class :character
## Mode :character
##
##
## HW2_cit$age_sq=HW2_cit$age_centered^2 #Using nlme to fit No Change Model to CIT data #Model 1: Unconditional Means
cit_no_change_A <- nlme(cit ~ b_1i,
data = HW2_cit,
fixed = b_1i~1,
random = b_1i~1,
groups=~ID,
start=c(18),
method="ML",
na.action=na.omit)
summary(cit_no_change_A)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: cit ~ b_1i
## Data: HW2_cit
## AIC BIC logLik
## 211.1702 225.224 -102.5851
##
## Random effects:
## Formula: b_1i ~ 1 | ID
## b_1i Residual
## StdDev: 0.0369354 0.2726513
##
## Fixed effects: b_1i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 3.001343 0.009993465 600 300.3305 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.16953741 -0.67275777 -0.02432779 0.68622615 2.96067594
##
## Number of Observations: 800
## Number of Groups: 200#Model 2: Unconditional Linear Growth
cit_linear <- nlme(cit ~ b_1i+b_2i*((time)/1),
data = HW2_cit,
fixed=b_1i+b_2i~1,
random=b_1i+b_2i~1,
groups=~ID,
start=c(18,.04),
method="ML",
na.action=na.omit)## Warning in nlme.formula(cit ~ b_1i + b_2i * ((time)/1), data = HW2_cit, :
## Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase
## 'msMaxIter'!summary(cit_linear)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: cit ~ b_1i + b_2i * ((time)/1)
## Data: HW2_cit
## AIC BIC logLik
## 211.4721 239.5798 -99.73605
##
## Random effects:
## Formula: list(b_1i ~ 1, b_2i ~ 1)
## Level: ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b_1i 0.05359027 b_1i
## b_2i 0.04389842 -0.999
## Residual 0.26480764
##
## Fixed effects: b_1i + b_2i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 2.9937884 0.023273112 599 128.63722 0.0000
## b_2i 0.0030217 0.008941953 599 0.33792 0.7355
## Correlation:
## b_1i
## b_2i -0.901
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.9645528 -0.6562038 -0.0468912 0.6840902 2.9782168
##
## Number of Observations: 800
## Number of Groups: 200#Model 3: Unconditional Quadratic
cit_quadratic <- nlme(cit ~ b_1i+b_2i*(time/1)+b_3i*(time/1)^2,
data = HW2_cit,
fixed=b_1i+b_2i+b_3i~1,
random=b_1i+b_2i+b_3i~1,
groups=~ID,
start=c(18,.04,1),
method="ML",
na.action=na.omit)## Warning in nlme.formula(cit ~ b_1i + b_2i * (time/1) + b_3i * (time/1)^2, :
## Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase
## 'msMaxIter'!summary(cit_quadratic)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: cit ~ b_1i + b_2i * (time/1) + b_3i * (time/1)^2
## Data: HW2_cit
## AIC BIC logLik
## 215.6894 262.5355 -97.8447
##
## Random effects:
## Formula: list(b_1i ~ 1, b_2i ~ 1, b_3i ~ 1)
## Level: ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## b_1i 2.077931e-05 b_1i b_2i
## b_2i 8.552288e-03 0.929
## b_3i 1.035532e-02 -0.963 -0.973
## Residual 2.638964e-01
##
## Fixed effects: b_1i + b_2i + b_3i ~ 1
## Value Std.Error DF t-value p-value
## b_1i 2.9042129 0.05204573 598 55.80118 0.0000
## b_2i 0.0926019 0.04748424 598 1.95016 0.0516
## b_3i -0.0179167 0.00937643 598 -1.91082 0.0565
## Correlation:
## b_1i b_2i
## b_2i -0.955
## b_3i 0.895 -0.982
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.11101024 -0.64250606 -0.05673723 0.67537729 2.96505275
##
## Number of Observations: 800
## Number of Groups: 200rand_intercept1 <- lmer(cit ~ 1 + time + age_centered + + time*age_centered + time*
+ (1 | ID),
data = HW2_long,
REML = FALSE)
summary(rand_intercept1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: cit ~ 1 + time + age_centered + +time * age_centered + time *
## +(1 | ID)
## Data: HW2_long
##
## AIC BIC logLik deviance df.resid
## 214.1 242.2 -101.1 202.1 794
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2784 -0.6735 -0.0283 0.6835 2.9394
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.001334 0.03653
## Residual 0.074078 0.27217
## Number of obs: 800, groups: ID, 200
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.030435 0.048398 62.615
## time -0.018584 0.017567 -1.058
## age_centered -0.002576 0.002966 -0.868
## time:age_centered 0.001519 0.001077 1.411
##
## Correlation of Fixed Effects:
## (Intr) time ag_cnt
## time -0.907
## age_centerd -0.872 0.791
## tim:g_cntrd 0.791 -0.872 -0.907# Random slope 1
rand_slope1 <- lmer(cit ~ 1 + time + age_centered + + time*age_centered +
+ (1 + time | ID),
data = HW2_long,
REML = FALSE)## boundary (singular) fit: see ?isSingularsummary(rand_slope1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: cit ~ 1 + time + age_centered + +time * age_centered + +(1 +
## time | ID)
## Data: HW2_long
##
## AIC BIC logLik deviance df.resid
## 213.0 250.5 -98.5 197.0 792
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.05480 -0.65223 -0.06621 0.68987 2.96980
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.002580 0.05080
## time 0.001801 0.04243 -1.00
## Residual 0.070105 0.26477
## Number of obs: 800, groups: ID, 200
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.030435 0.047373 63.970
## time -0.018584 0.018154 -1.024
## age_centered -0.002576 0.002903 -0.887
## time:age_centered 0.001519 0.001113 1.365
##
## Correlation of Fixed Effects:
## (Intr) time ag_cnt
## time -0.901
## age_centerd -0.872 0.786
## tim:g_cntrd 0.786 -0.872 -0.901
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularrand_intercept2 <- lmer(cit ~ 1 + time + age_centered + Relationship_status + time*age_centered + time*Relationship_status
+ (1 | ID),
data = HW2_long,
REML = FALSE)
summary(rand_intercept2)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: cit ~ 1 + time + age_centered + Relationship_status + time *
## age_centered + time * Relationship_status + (1 | ID)
## Data: HW2_long
##
## AIC BIC logLik deviance df.resid
## 216.6 254.1 -100.3 200.6 792
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2656 -0.6604 -0.0274 0.7066 2.9381
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.001367 0.03698
## Residual 0.073905 0.27185
## Number of obs: 800, groups: ID, 200
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.997501 0.059667 50.237
## time -0.003530 0.021654 -0.163
## age_centered -0.002746 0.002968 -0.925
## Relationship_status 0.009660 0.010256 0.942
## time:age_centered 0.001596 0.001077 1.482
## time:Relationship_status -0.004416 0.003722 -1.186
##
## Correlation of Fixed Effects:
## (Intr) time ag_cnt Rltns_ tm:g_c
## time -0.907
## age_centerd -0.669 0.607
## Rltnshp_stt -0.586 0.532 -0.061
## tim:g_cntrd 0.607 -0.669 -0.907 0.055
## tm:Rltnshp_ 0.532 -0.586 0.055 -0.907 -0.061# Random slope 2
rand_slope2 <- lmer(cit ~ 1 + time + age_centered + Relationship_status + time*age_centered + time*Relationship_status
+ (1 + time | ID),
data = HW2_long,
REML = FALSE)## boundary (singular) fit: see ?isSingularsummary(rand_slope2)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: cit ~ 1 + time + age_centered + Relationship_status + time *
## age_centered + time * Relationship_status + (1 + time | ID)
## Data: HW2_long
##
## AIC BIC logLik deviance df.resid
## 215.7 262.5 -97.8 195.7 790
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0483 -0.6519 -0.0637 0.6912 2.9649
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.002386 0.04885
## time 0.001735 0.04165 -1.00
## Residual 0.070046 0.26466
## Number of obs: 800, groups: ID, 200
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.997501 0.058386 51.339
## time -0.003530 0.022349 -0.158
## age_centered -0.002746 0.002905 -0.945
## Relationship_status 0.009660 0.010036 0.963
## time:age_centered 0.001596 0.001112 1.436
## time:Relationship_status -0.004416 0.003841 -1.149
##
## Correlation of Fixed Effects:
## (Intr) time ag_cnt Rltns_ tm:g_c
## time -0.901
## age_centerd -0.669 0.603
## Rltnshp_stt -0.586 0.528 -0.061
## tim:g_cntrd 0.603 -0.669 -0.901 0.055
## tm:Rltnshp_ 0.528 -0.586 0.055 -0.901 -0.061
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular