library(dplyr)
## Warning: package 'dplyr' was built under R version 3.6.3
##
## 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
library(car) # leveneTest, vif
## Warning: package 'car' was built under R version 3.6.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.6.1
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
library(MASS) # lda
## Warning: package 'MASS' was built under R version 3.6.3
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library(psy) # cronbach
library(latexpdf)
library(mediation)
## Warning: package 'mediation' was built under R version 3.6.3
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 3.6.3
## Loading required package: mvtnorm
## Warning: package 'mvtnorm' was built under R version 3.6.2
## Loading required package: sandwich
## Warning: package 'sandwich' was built under R version 3.6.1
## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car
## mediation: Causal Mediation Analysis
## Version: 4.5.0
library(knitr)
## Warning: package 'knitr' was built under R version 3.6.3
library(boot)
##
## Attaching package: 'boot'
## The following object is masked from 'package:car':
##
## logit
library(haven)
## Warning: package 'haven' was built under R version 3.6.3
library(bda)
## Warning: package 'bda' was built under R version 3.6.3
## bda v14 (Bin Wang, 2020)
library(reshape2)
## Warning: package 'reshape2' was built under R version 3.6.3
library(Hmisc) # correlation matrix
## Warning: package 'Hmisc' was built under R version 3.6.3
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 3.6.3
##
## Attaching package: 'lattice'
## The following object is masked from 'package:boot':
##
## melanoma
## Loading required package: survival
## Warning: package 'survival' was built under R version 3.6.3
##
## Attaching package: 'survival'
## The following object is masked from 'package:boot':
##
## aml
## Loading required package: Formula
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.6.3
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
##
## src, summarize
## The following objects are masked from 'package:base':
##
## format.pval, units
library(e1071)
## Warning: package 'e1071' was built under R version 3.6.3
##
## Attaching package: 'e1071'
## The following object is masked from 'package:Hmisc':
##
## impute
##################
## STUDY 1 ####
##################
cat("\014") # cleans screen
rm(list=ls(all=T)) # remove variables in working memory
setwd("C:/Users/ev193805/Downloads/Papers Pipeline 2/In-Progress/Marketing Letters/Reto")
study1=read_sav("Study 1 V2 2020-06-23.sav")
Item1AH<-melt(study1,id.vars=c("ID"),measure.vars=c("Q1_1","Q7_1","Q13_1"),variable.name="Appearance Homophily", value.name="Item1AH")
Item2AH<-melt(study1,id.vars=c("ID"),measure.vars=c("Q1_2","Q7_2","Q13_2"),variable.name="Appearance Homophily", value.name="Item2AH")
Item3AH<-melt(study1,id.vars=c("ID"),measure.vars=c("Q1_3","Q7_3","Q13_3"),variable.name="Appearance Homophily", value.name="Item3AH")
Item4AH<-melt(study1,id.vars=c("ID"),measure.vars=c("Q1_4","Q7_4","Q13_4"),variable.name="Appearance Homophily", value.name="Item4AH")
Item1Orig<-melt(study1,id.vars=c("ID"),measure.vars=c("Q3_1","Q9_1","Q15_1"),variable.name="Originality", value.name="Item1Orig")
Item2Orig<-melt(study1,id.vars=c("ID"),measure.vars=c("Q3_2","Q9_2","Q15_2"),variable.name="Originality", value.name="Item2Orig")
Item3Orig<-melt(study1,id.vars=c("ID"),measure.vars=c("Q3_3","Q9_3","Q15_3"),variable.name="Originality", value.name="Item3Orig")
Item4Orig<-melt(study1,id.vars=c("ID"),measure.vars=c("Q3_4","Q9_4","Q15_4"),variable.name="Originality", value.name="Item4Orig")
Item5Orig<-melt(study1,id.vars=c("ID"),measure.vars=c("Q3_5","Q9_5","Q15_5"),variable.name="Originality", value.name="Item5Orig")
Item1Nat<-melt(study1,id.vars=c("ID"),measure.vars=c("Q4_1","Q10_1","Q16_1"),variable.name="Naturalness", value.name="Item1Nat")
Item2Nat<-melt(study1,id.vars=c("ID"),measure.vars=c("Q4_2","Q10_2","Q16_2"),variable.name="Naturalness", value.name="Item2Nat")
Item3Nat<-melt(study1,id.vars=c("ID"),measure.vars=c("Q4_3","Q10_3","Q16_3"),variable.name="Naturalness", value.name="Item3Nat")
Item1Rel<-melt(study1,id.vars=c("ID"),measure.vars=c("Q31_1","Q32_1","Q33_1"),variable.name="Reliability", value.name="Item1Rel")
Item2Rel<-melt(study1,id.vars=c("ID"),measure.vars=c("Q31_2","Q32_2","Q33_2"),variable.name="Reliability", value.name="Item2Rel")
Item3Rel<-melt(study1,id.vars=c("ID"),measure.vars=c("Q31_3","Q32_3","Q33_3"),variable.name="Reliability", value.name="Item3Rel")
Item4Rel<-melt(study1,id.vars=c("ID"),measure.vars=c("Q31_4","Q32_4","Q33_4"),variable.name="Reliability", value.name="Item4Rel")
Item5Rel<-melt(study1,id.vars=c("ID"),measure.vars=c("Q31_5","Q32_5","Q33_5"),variable.name="Reliability", value.name="Item5Rel")
Item6Rel<-melt(study1,id.vars=c("ID"),measure.vars=c("Q31_6","Q32_6","Q33_6"),variable.name="Reliability", value.name="Item6Rel")
Item7Rel<-melt(study1,id.vars=c("ID"),measure.vars=c("Q31_7","Q32_7","Q33_7"),variable.name="Reliability", value.name="Item7Rel")
Item8Rel<-melt(study1,id.vars=c("ID"),measure.vars=c("Q31_8","Q32_8","Q33_8"),variable.name="Reliability", value.name="Item8Rel")
Item9Rel<-melt(study1,id.vars=c("ID"),measure.vars=c("Q31_9","Q32_9","Q33_9"),variable.name="Reliability", value.name="Item9Rel")
study1melt=cbind(Item1AH[1],
Item1AH[3],Item2AH[3],Item3AH[3],Item4AH[3],
Item1Orig[3],Item2Orig[3],Item3Orig[3],Item4Orig[3],Item5Orig[3],
Item1Nat[3],Item2Nat[3],Item3Nat[3],
Item1Rel[3],Item2Rel[3],Item3Rel[3],Item4Rel[3],Item5Rel[3],
Item6Rel[3],Item7Rel[3],Item8Rel[3],Item9Rel[3])
summary(study1melt)
## ID Item1AH Item2AH Item3AH Item4AH
## Min. :248.0 Min. :1.0 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:340.0 1st Qu.:3.0 1st Qu.:1.000 1st Qu.:3.000 1st Qu.:1.000
## Median :437.0 Median :5.0 Median :3.000 Median :5.000 Median :3.000
## Mean :477.3 Mean :4.7 Mean :3.327 Mean :4.633 Mean :3.286
## 3rd Qu.:625.0 3rd Qu.:7.0 3rd Qu.:5.000 3rd Qu.:6.000 3rd Qu.:5.000
## Max. :792.0 Max. :7.0 Max. :7.000 Max. :7.000 Max. :7.000
## Item1Orig Item2Orig Item3Orig Item4Orig Item5Orig
## Min. :1.000 Min. :1.000 Min. :1.00 Min. :1.000 Min. :1.000
## 1st Qu.:3.000 1st Qu.:3.000 1st Qu.:2.00 1st Qu.:2.000 1st Qu.:3.000
## Median :4.000 Median :4.000 Median :4.00 Median :4.000 Median :4.000
## Mean :4.173 Mean :3.856 Mean :3.66 Mean :3.796 Mean :3.908
## 3rd Qu.:5.000 3rd Qu.:5.000 3rd Qu.:5.00 3rd Qu.:5.000 3rd Qu.:5.000
## Max. :7.000 Max. :7.000 Max. :7.00 Max. :7.000 Max. :7.000
## Item1Nat Item2Nat Item3Nat Item1Rel
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:3.000 1st Qu.:3.000 1st Qu.:2.000 1st Qu.:4.000
## Median :5.000 Median :4.000 Median :4.000 Median :4.000
## Mean :4.541 Mean :4.297 Mean :3.751 Mean :4.506
## 3rd Qu.:6.000 3rd Qu.:6.000 3rd Qu.:5.000 3rd Qu.:5.000
## Max. :7.000 Max. :7.000 Max. :7.000 Max. :7.000
## Item2Rel Item3Rel Item4Rel Item5Rel
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:4.000 1st Qu.:4.000 1st Qu.:4.000 1st Qu.:4.000
## Median :4.000 Median :4.000 Median :4.000 Median :4.000
## Mean :4.515 Mean :4.539 Mean :4.512 Mean :4.519
## 3rd Qu.:5.000 3rd Qu.:5.000 3rd Qu.:5.000 3rd Qu.:5.000
## Max. :7.000 Max. :7.000 Max. :7.000 Max. :7.000
## Item6Rel Item7Rel Item8Rel Item9Rel
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:4.000 1st Qu.:4.000 1st Qu.:4.000 1st Qu.:4.000
## Median :4.000 Median :4.000 Median :4.000 Median :5.000
## Mean :4.532 Mean :4.465 Mean :4.596 Mean :4.677
## 3rd Qu.:6.000 3rd Qu.:5.000 3rd Qu.:6.000 3rd Qu.:6.000
## Max. :7.000 Max. :7.000 Max. :7.000 Max. :7.000
rcorr(as.matrix(study1melt[,2:22]))
## Item1AH Item2AH Item3AH Item4AH Item1Orig Item2Orig Item3Orig
## Item1AH 1.00 -0.36 0.72 -0.44 0.03 -0.01 -0.04
## Item2AH -0.36 1.00 -0.33 0.65 0.19 0.16 0.20
## Item3AH 0.72 -0.33 1.00 -0.38 -0.01 -0.04 -0.04
## Item4AH -0.44 0.65 -0.38 1.00 0.22 0.20 0.22
## Item1Orig 0.03 0.19 -0.01 0.22 1.00 0.62 0.50
## Item2Orig -0.01 0.16 -0.04 0.20 0.62 1.00 0.64
## Item3Orig -0.04 0.20 -0.04 0.22 0.50 0.64 1.00
## Item4Orig -0.02 0.20 -0.03 0.21 0.60 0.72 0.64
## Item5Orig 0.01 0.17 0.00 0.20 0.66 0.69 0.61
## Item1Nat 0.03 0.15 0.03 0.15 0.42 0.34 0.16
## Item2Nat 0.00 0.19 0.03 0.20 0.40 0.33 0.20
## Item3Nat 0.09 0.00 0.09 0.01 -0.17 -0.06 0.05
## Item1Rel -0.01 0.18 0.00 0.22 0.38 0.24 0.13
## Item2Rel -0.03 0.19 -0.02 0.23 0.39 0.26 0.16
## Item3Rel 0.00 0.20 0.02 0.24 0.44 0.28 0.17
## Item4Rel -0.03 0.20 -0.04 0.26 0.42 0.26 0.18
## Item5Rel 0.00 0.20 -0.01 0.27 0.40 0.27 0.18
## Item6Rel -0.03 0.23 -0.02 0.25 0.40 0.27 0.19
## Item7Rel 0.00 0.18 0.01 0.24 0.39 0.27 0.14
## Item8Rel 0.00 0.19 0.02 0.21 0.38 0.26 0.19
## Item9Rel -0.02 0.17 0.02 0.19 0.35 0.20 0.17
## Item4Orig Item5Orig Item1Nat Item2Nat Item3Nat Item1Rel Item2Rel
## Item1AH -0.02 0.01 0.03 0.00 0.09 -0.01 -0.03
## Item2AH 0.20 0.17 0.15 0.19 0.00 0.18 0.19
## Item3AH -0.03 0.00 0.03 0.03 0.09 0.00 -0.02
## Item4AH 0.21 0.20 0.15 0.20 0.01 0.22 0.23
## Item1Orig 0.60 0.66 0.42 0.40 -0.17 0.38 0.39
## Item2Orig 0.72 0.69 0.34 0.33 -0.06 0.24 0.26
## Item3Orig 0.64 0.61 0.16 0.20 0.05 0.13 0.16
## Item4Orig 1.00 0.69 0.30 0.28 -0.02 0.21 0.22
## Item5Orig 0.69 1.00 0.34 0.33 -0.05 0.29 0.31
## Item1Nat 0.30 0.34 1.00 0.69 -0.41 0.47 0.47
## Item2Nat 0.28 0.33 0.69 1.00 -0.28 0.45 0.42
## Item3Nat -0.02 -0.05 -0.41 -0.28 1.00 -0.20 -0.19
## Item1Rel 0.21 0.29 0.47 0.45 -0.20 1.00 0.87
## Item2Rel 0.22 0.31 0.47 0.42 -0.19 0.87 1.00
## Item3Rel 0.25 0.34 0.48 0.43 -0.21 0.83 0.83
## Item4Rel 0.26 0.32 0.44 0.40 -0.16 0.76 0.81
## Item5Rel 0.24 0.32 0.45 0.40 -0.14 0.77 0.80
## Item6Rel 0.25 0.30 0.46 0.41 -0.18 0.77 0.80
## Item7Rel 0.24 0.27 0.47 0.40 -0.16 0.74 0.76
## Item8Rel 0.25 0.30 0.43 0.36 -0.16 0.72 0.74
## Item9Rel 0.23 0.26 0.44 0.36 -0.19 0.68 0.70
## Item3Rel Item4Rel Item5Rel Item6Rel Item7Rel Item8Rel Item9Rel
## Item1AH 0.00 -0.03 0.00 -0.03 0.00 0.00 -0.02
## Item2AH 0.20 0.20 0.20 0.23 0.18 0.19 0.17
## Item3AH 0.02 -0.04 -0.01 -0.02 0.01 0.02 0.02
## Item4AH 0.24 0.26 0.27 0.25 0.24 0.21 0.19
## Item1Orig 0.44 0.42 0.40 0.40 0.39 0.38 0.35
## Item2Orig 0.28 0.26 0.27 0.27 0.27 0.26 0.20
## Item3Orig 0.17 0.18 0.18 0.19 0.14 0.19 0.17
## Item4Orig 0.25 0.26 0.24 0.25 0.24 0.25 0.23
## Item5Orig 0.34 0.32 0.32 0.30 0.27 0.30 0.26
## Item1Nat 0.48 0.44 0.45 0.46 0.47 0.43 0.44
## Item2Nat 0.43 0.40 0.40 0.41 0.40 0.36 0.36
## Item3Nat -0.21 -0.16 -0.14 -0.18 -0.16 -0.16 -0.19
## Item1Rel 0.83 0.76 0.77 0.77 0.74 0.72 0.68
## Item2Rel 0.83 0.81 0.80 0.80 0.76 0.74 0.70
## Item3Rel 1.00 0.82 0.81 0.77 0.75 0.74 0.70
## Item4Rel 0.82 1.00 0.86 0.81 0.79 0.75 0.73
## Item5Rel 0.81 0.86 1.00 0.81 0.79 0.78 0.71
## Item6Rel 0.77 0.81 0.81 1.00 0.81 0.78 0.76
## Item7Rel 0.75 0.79 0.79 0.81 1.00 0.77 0.70
## Item8Rel 0.74 0.75 0.78 0.78 0.77 1.00 0.81
## Item9Rel 0.70 0.73 0.71 0.76 0.70 0.81 1.00
##
## n= 963
##
##
## P
## Item1AH Item2AH Item3AH Item4AH Item1Orig Item2Orig Item3Orig
## Item1AH 0.0000 0.0000 0.0000 0.3951 0.6986 0.1735
## Item2AH 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item3AH 0.0000 0.0000 0.0000 0.7410 0.1876 0.2562
## Item4AH 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item1Orig 0.3951 0.0000 0.7410 0.0000 0.0000 0.0000
## Item2Orig 0.6986 0.0000 0.1876 0.0000 0.0000 0.0000
## Item3Orig 0.1735 0.0000 0.2562 0.0000 0.0000 0.0000
## Item4Orig 0.4543 0.0000 0.3723 0.0000 0.0000 0.0000 0.0000
## Item5Orig 0.8245 0.0000 0.8979 0.0000 0.0000 0.0000 0.0000
## Item1Nat 0.3785 0.0000 0.3815 0.0000 0.0000 0.0000 0.0000
## Item2Nat 0.9452 0.0000 0.3831 0.0000 0.0000 0.0000 0.0000
## Item3Nat 0.0062 0.9380 0.0046 0.8404 0.0000 0.0559 0.1614
## Item1Rel 0.6517 0.0000 0.9553 0.0000 0.0000 0.0000 0.0000
## Item2Rel 0.3478 0.0000 0.5112 0.0000 0.0000 0.0000 0.0000
## Item3Rel 0.9778 0.0000 0.6020 0.0000 0.0000 0.0000 0.0000
## Item4Rel 0.2829 0.0000 0.2150 0.0000 0.0000 0.0000 0.0000
## Item5Rel 0.9809 0.0000 0.6819 0.0000 0.0000 0.0000 0.0000
## Item6Rel 0.3867 0.0000 0.5110 0.0000 0.0000 0.0000 0.0000
## Item7Rel 0.9563 0.0000 0.6761 0.0000 0.0000 0.0000 0.0000
## Item8Rel 0.9219 0.0000 0.5990 0.0000 0.0000 0.0000 0.0000
## Item9Rel 0.4670 0.0000 0.6390 0.0000 0.0000 0.0000 0.0000
## Item4Orig Item5Orig Item1Nat Item2Nat Item3Nat Item1Rel Item2Rel
## Item1AH 0.4543 0.8245 0.3785 0.9452 0.0062 0.6517 0.3478
## Item2AH 0.0000 0.0000 0.0000 0.0000 0.9380 0.0000 0.0000
## Item3AH 0.3723 0.8979 0.3815 0.3831 0.0046 0.9553 0.5112
## Item4AH 0.0000 0.0000 0.0000 0.0000 0.8404 0.0000 0.0000
## Item1Orig 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item2Orig 0.0000 0.0000 0.0000 0.0000 0.0559 0.0000 0.0000
## Item3Orig 0.0000 0.0000 0.0000 0.0000 0.1614 0.0000 0.0000
## Item4Orig 0.0000 0.0000 0.0000 0.6337 0.0000 0.0000
## Item5Orig 0.0000 0.0000 0.0000 0.1014 0.0000 0.0000
## Item1Nat 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item2Nat 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item3Nat 0.6337 0.1014 0.0000 0.0000 0.0000 0.0000
## Item1Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item2Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item3Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item4Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item5Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item6Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item7Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item8Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item9Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item3Rel Item4Rel Item5Rel Item6Rel Item7Rel Item8Rel Item9Rel
## Item1AH 0.9778 0.2829 0.9809 0.3867 0.9563 0.9219 0.4670
## Item2AH 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item3AH 0.6020 0.2150 0.6819 0.5110 0.6761 0.5990 0.6390
## Item4AH 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item1Orig 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item2Orig 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item3Orig 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item4Orig 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item5Orig 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item1Nat 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item2Nat 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item3Nat 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item1Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item2Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item3Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item4Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item5Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item6Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item7Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item8Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## Item9Rel 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
cronbach(cbind(Item1AH[3],-Item2AH[3],Item3AH[3],-Item4AH[3])) ## .79
## $sample.size
## [1] 963
##
## $number.of.items
## [1] 4
##
## $alpha
## [1] 0.7874914
cronbach(cbind(Item1Orig[3],Item2Orig[3],Item3Orig[3],Item4Orig[3],Item5Orig[3])) ## .90
## $sample.size
## [1] 963
##
## $number.of.items
## [1] 5
##
## $alpha
## [1] 0.8982091
cronbach(cbind(Item1Nat[3],Item2Nat[3],-Item3Nat[3])) ## .72
## $sample.size
## [1] 963
##
## $number.of.items
## [1] 3
##
## $alpha
## [1] 0.715038
cronbach(cbind(Item1Rel[3],Item2Rel[3],Item3Rel[3],Item4Rel[3],Item5Rel[3],
Item6Rel[3],Item7Rel[3],Item8Rel[3],Item9Rel[3])) ## .97
## $sample.size
## [1] 963
##
## $number.of.items
## [1] 9
##
## $alpha
## [1] 0.9681528
study1meltwithoutrev=cbind(Item1AH[1],
Item1AH[3],-Item2AH[3],Item3AH[3],-Item4AH[3],
Item1Orig[3],Item2Orig[3],Item3Orig[3],Item4Orig[3],Item5Orig[3],
Item1Nat[3],Item2Nat[3],-Item3Nat[3],
Item1Rel[3],Item2Rel[3],Item3Rel[3],Item4Rel[3],Item5Rel[3],
Item6Rel[3],Item7Rel[3],Item8Rel[3],Item9Rel[3])
factanal(study1meltwithoutrev[2:22],4,rotation="varimax") ## four components explain 67% of the variance
##
## Call:
## factanal(x = study1meltwithoutrev[2:22], factors = 4, rotation = "varimax")
##
## Uniquenesses:
## Item1AH Item2AH Item3AH Item4AH Item1Orig Item2Orig Item3Orig Item4Orig
## 0.241 0.711 0.350 0.615 0.403 0.290 0.411 0.290
## Item5Orig Item1Nat Item2Nat Item3Nat Item1Rel Item2Rel Item3Rel Item4Rel
## 0.308 0.161 0.433 0.772 0.226 0.186 0.197 0.174
## Item5Rel Item6Rel Item7Rel Item8Rel Item9Rel
## 0.179 0.201 0.249 0.270 0.343
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4
## Item1AH 0.869
## Item2AH -0.153 -0.163 0.483
## Item3AH 0.804
## Item4AH -0.196 -0.183 0.556
## Item1Orig 0.310 0.677 0.203
## Item2Orig 0.134 0.820 0.123
## Item3Orig 0.757
## Item4Orig 0.121 0.827
## Item5Orig 0.196 0.800 0.112
## Item1Nat 0.343 0.235 0.816
## Item2Nat 0.310 0.255 0.637
## Item3Nat 0.118 0.455
## Item1Rel 0.846 0.113 0.204
## Item2Rel 0.871 0.127 0.181
## Item3Rel 0.859 0.165 0.187
## Item4Rel 0.881 0.160 0.118
## Item5Rel 0.880 0.156 0.132
## Item6Rel 0.862 0.151 0.158
## Item7Rel 0.837 0.137 0.173
## Item8Rel 0.829 0.160 0.129
## Item9Rel 0.783 0.125 0.160
##
## Factor1 Factor2 Factor3 Factor4
## SS loadings 6.975 3.401 2.012 1.602
## Proportion Var 0.332 0.162 0.096 0.076
## Cumulative Var 0.332 0.494 0.590 0.666
##
## Test of the hypothesis that 4 factors are sufficient.
## The chi square statistic is 1108.59 on 132 degrees of freedom.
## The p-value is 5.66e-154
summary(prcomp(study1meltwithoutrev[2:22])) ## four components explain 72% of the variance
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 4.4103 2.9453 2.726 2.00797 1.78000 1.37267 1.15327
## Proportion of Variance 0.3533 0.1576 0.135 0.07324 0.05756 0.03423 0.02416
## Cumulative Proportion 0.3533 0.5109 0.646 0.71919 0.77675 0.81098 0.83514
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 1.14034 1.0255 1.00179 0.96586 0.92047 0.88396 0.84372
## Proportion of Variance 0.02362 0.0191 0.01823 0.01695 0.01539 0.01419 0.01293
## Cumulative Proportion 0.85876 0.8779 0.89610 0.91304 0.92844 0.94263 0.95556
## PC15 PC16 PC17 PC18 PC19 PC20 PC21
## Standard deviation 0.74598 0.66694 0.62180 0.57015 0.52733 0.49581 0.4575
## Proportion of Variance 0.01011 0.00808 0.00702 0.00591 0.00505 0.00447 0.0038
## Cumulative Proportion 0.96567 0.97375 0.98078 0.98668 0.99173 0.99620 1.0000
screeplot(prcomp(study1meltwithoutrev[2:22]),type="lines")
biplot(prcomp(study1meltwithoutrev[2:22],scale.=T),cex=0.5,xlabs=rep(".",nrow(study1meltwithoutrev[2:22])))
## BMI Difference
BMIDiff<-melt(study1,id.vars=c("ID"),measure.vars=c("B_17_ABS","B_22_ABS","B_28_ABS"),variable.name="BMI Difference", value.name="Value")
## Warning: attributes are not identical across measure variables; they will be
## dropped
study1$RaceDiff=ifelse(study1$Q31==3,0,1)
RaceDiff<-melt(study1,id.vars=c("ID"),measure.vars=c("RaceDiff","RaceDiff","RaceDiff"),variable.name="Race Match", value.name="Dummy")
study1$GenDiff=ifelse(study1$Q27==1,0,1)
GenDiff<-melt(study1,id.vars=c("ID"),measure.vars=c("GenDiff","GenDiff","GenDiff"),variable.name="Gen Match", value.name="Dummy")
Age<-melt(study1,id.vars=c("ID"),measure.vars=c("Q26","Q26","Q26"),variable.name="Age", value.name="Year of birth")
study1meltwithoutrev$BMIDiff=BMIDiff$Value
study1meltwithoutrev$RaceDiff=RaceDiff$Dummy
study1meltwithoutrev$GenDiff=GenDiff$Dummy
study1meltwithoutrev$AgeDiff=abs(22-2019-Age$`Year of birth`)
## Mixed Mixed Effects Modeling
## https://ademos.people.uic.edu/Chapter17.html#11_example:_national_pizza_study
library (lmerTest) # Mixed model package by Douglas Bates, comes w/ pvalues!
## Warning: package 'lmerTest' was built under R version 3.6.3
## Loading required package: lme4
## Warning: package 'lme4' was built under R version 3.6.3
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library (texreg) #Helps us make tables of the mixed models
## Warning: package 'texreg' was built under R version 3.6.3
## Version: 1.37.5
## Date: 2020-06-17
## Author: Philip Leifeld (University of Essex)
##
## Consider submitting praise using the praise or praise_interactive functions.
## Please cite the JSS article in your publications -- see citation("texreg").
library (afex) # Easy ANOVA package to compare model fits
## Warning: package 'afex' was built under R version 3.6.3
## ************
## Welcome to afex. For support visit: http://afex.singmann.science/
## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'KR', 'S', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - NEWS: library('emmeans') now needs to be called explicitly!
## - Get and set global package options with: afex_options()
## - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************
##
## Attaching package: 'afex'
## The following object is masked from 'package:lme4':
##
## lmer
library (plyr) # Data manipulator package
## Warning: package 'plyr' was built under R version 3.6.3
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:Hmisc':
##
## is.discrete, summarize
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
library (ggplot2) # GGplot package for visualizing data
study1meltwithoutrev$AdType=0
study1meltwithoutrev[c(1:321),27]=1
study1meltwithoutrev[c(322:642),27]=2
study1meltwithoutrev[c(643:963),27]=3
study1meltwithoutrev$AvgAH=(study1meltwithoutrev$Item1AH-study1meltwithoutrev$Item2AH
+study1meltwithoutrev$Item3AH-study1meltwithoutrev$Item4AH)/4
study1meltwithoutrev$AvgOrig=(study1meltwithoutrev$Item1Orig+study1meltwithoutrev$Item2Orig
+study1meltwithoutrev$Item3Orig+study1meltwithoutrev$Item4Orig+study1meltwithoutrev$Item5Orig)/5
## Only Caucasian women
study1b=subset(study1meltwithoutrev,study1meltwithoutrev$GenDiff==0)
study1b=subset(study1b,study1b$RaceDiff==0)
## Model fitness
nullmodel1<-lmer(AvgAH~1+(1|ID),study1b,REML=FALSE)
nullmodel2<-lmer(AvgAH~1+(1+BMIDiff|ID),study1b,REML=FALSE)
## boundary (singular) fit: see ?isSingular
nullmodel3<-lmer(AvgAH~1+(1+AgeDiff|ID),study1b,REML=FALSE)
## boundary (singular) fit: see ?isSingular
nullmodel4<-lmer(AvgAH~1+(1+GenDiff|ID),study1b,REML=FALSE)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Hessian is numerically singular: parameters are not uniquely determined
## Warning in as_lmerModLT(model, devfun): Model may not have converged with 2
## eigenvalues close to zero: -1.1e-30 -2.7e-27
nullmodel5<-lmer(AvgAH~1+(1+AdType|ID),study1b,REML=FALSE)
anova(nullmodel1,nullmodel2,nullmodel3,nullmodel4,nullmodel5)
## Data: study1b
## Models:
## nullmodel1: AvgAH ~ 1 + (1 | ID)
## nullmodel2: AvgAH ~ 1 + (1 + BMIDiff | ID)
## nullmodel3: AvgAH ~ 1 + (1 + AgeDiff | ID)
## nullmodel4: AvgAH ~ 1 + (1 + GenDiff | ID)
## nullmodel5: AvgAH ~ 1 + (1 + AdType | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nullmodel1 3 1029.2 1041.6 -511.63 1023.2
## nullmodel2 5 1025.8 1046.3 -507.91 1015.8 7.4338 2 0.02431 *
## nullmodel3 5 1033.2 1053.7 -511.61 1023.2 0.0000 0 1.00000
## nullmodel4 5 1033.2 1053.8 -511.63 1023.2 0.0000 0 1.00000
## nullmodel5 5 1032.1 1052.7 -511.07 1022.1 1.1076 0 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Fixed effects
m1=lmer(AvgAH~AdType+(1+AdType|ID),study1b,REML=FALSE)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00349114 (tol = 0.002, component 1)
summary(m1)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: AvgAH ~ AdType + (1 + AdType | ID)
## Data: study1b
##
## AIC BIC logLik deviance df.resid
## 1033.7 1058.3 -510.8 1021.7 441
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0351 -0.2843 0.0570 0.3943 2.9456
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.19975 0.4469
## AdType 0.01849 0.1360 -0.27
## Residual 0.41214 0.6420
## Number of obs: 447, groups: ID, 149
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.95973 0.08829 149.13664 44.850 <2e-16 ***
## AdType -0.02685 0.03882 149.15272 -0.692 0.49
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## AdType -0.839
## convergence code: 0
## Model failed to converge with max|grad| = 0.00349114 (tol = 0.002, component 1)
m2=lmer(AvgAH~BMIDiff+AdType+(1+AdType|ID),study1b,REML=FALSE)
summary(m2)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: AvgAH ~ BMIDiff + AdType + (1 + AdType | ID)
## Data: study1b
##
## AIC BIC logLik deviance df.resid
## 1034.7 1063.4 -510.3 1020.7 440
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0780 -0.3019 0.0441 0.4007 2.9189
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.18036 0.4247
## AdType 0.01472 0.1213 -0.14
## Residual 0.41292 0.6426
## Number of obs: 447, groups: ID, 149
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.043371 0.120759 204.871844 33.483 <2e-16 ***
## BMIDiff -0.006814 0.006770 255.651538 -1.006 0.315
## AdType -0.043391 0.041888 180.733405 -1.036 0.302
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) BMIDff
## BMIDiff -0.688
## AdType -0.827 0.392
m3=lmer(AvgAH~AgeDiff+BMIDiff+AdType+(1+AdType|ID),study1b,REML=FALSE)
summary(m3)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: AvgAH ~ AgeDiff + BMIDiff + AdType + (1 + AdType | ID)
## Data: study1b
##
## AIC BIC logLik deviance df.resid
## 1035.1 1067.9 -509.6 1019.1 439
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1192 -0.3109 0.0556 0.3962 2.8871
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.17292 0.4158
## AdType 0.01488 0.1220 -0.13
## Residual 0.41285 0.6425
## Number of obs: 447, groups: ID, 149
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -14.182534 14.489176 147.887119 -0.979 0.329
## AgeDiff 0.004584 0.003644 147.834856 1.258 0.210
## BMIDiff -0.006519 0.006745 255.496972 -0.967 0.335
## AdType -0.042675 0.041875 180.682238 -1.019 0.310
##
## Correlation of Fixed Effects:
## (Intr) AgeDff BMIDff
## AgeDiff -1.000
## BMIDiff -0.036 0.030
## AdType -0.019 0.012 0.391
anova(m1,m2,m3)
## Data: study1b
## Models:
## m1: AvgAH ~ AdType + (1 + AdType | ID)
## m2: AvgAH ~ BMIDiff + AdType + (1 + AdType | ID)
## m3: AvgAH ~ AgeDiff + BMIDiff + AdType + (1 + AdType | ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 6 1033.7 1058.3 -510.83 1021.7
## m2 7 1034.7 1063.4 -510.34 1020.7 0.9847 1 0.3211
## m3 8 1035.1 1067.9 -509.56 1019.1 1.5707 1 0.2101
## Regressions
study1Thin=subset(study1b,study1b$AdType==1)
study1Norm=subset(study1b,study1b$AdType==2)
study1Over=subset(study1b,study1b$AdType==3)
fit=lm(AvgAH~AgeDiff+BMIDiff,study1Thin)
summary(fit)
##
## Call:
## lm(formula = AvgAH ~ AgeDiff + BMIDiff, data = study1Thin)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9689 -0.1616 0.0941 0.3739 1.9064
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -28.901688 20.349878 -1.420 0.158
## AgeDiff 0.008281 0.005117 1.618 0.108
## BMIDiff -0.009055 0.009538 -0.949 0.344
##
## Residual standard error: 0.8185 on 146 degrees of freedom
## Multiple R-squared: 0.02547, Adjusted R-squared: 0.01212
## F-statistic: 1.908 on 2 and 146 DF, p-value: 0.1521
vif(fit)
## AgeDiff BMIDiff
## 1.007582 1.007582
shapiro.test(fit$res) ## Shapiro p > 0.01
##
## Shapiro-Wilk normality test
##
## data: fit$res
## W = 0.89164, p-value = 5.007e-09
hist(fit$res)
qqnorm(fit$res)
plot(fit$res)
skewness(fit$res) ## Skewness should be within the range +2 -2
## [1] -1.180623
kurtosis(fit$res) ## Kurtosis values should be within range of +7 -7
## [1] 2.597631
fit=lm(AvgAH~AgeDiff+BMIDiff,study1Norm)
summary(fit)
##
## Call:
## lm(formula = AvgAH ~ AgeDiff + BMIDiff, data = study1Norm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.92092 -0.19242 0.07227 0.33799 2.09342
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9033667 17.6193653 0.278 0.781
## AgeDiff -0.0002536 0.0044314 -0.057 0.954
## BMIDiff 0.0050620 0.0094613 0.535 0.593
##
## Residual standard error: 0.7109 on 146 degrees of freedom
## Multiple R-squared: 0.002, Adjusted R-squared: -0.01167
## F-statistic: 0.1463 on 2 and 146 DF, p-value: 0.864
vif(fit)
## AgeDiff BMIDiff
## 1.001758 1.001758
shapiro.test(fit$res) ## Shapiro p > 0.01
##
## Shapiro-Wilk normality test
##
## data: fit$res
## W = 0.91085, p-value = 6.176e-08
hist(fit$res)
qqnorm(fit$res)
plot(fit$res)
skewness(fit$res)
## [1] -0.8134035
kurtosis(fit$res)
## [1] 3.712247
fit=lm(AvgAH~AgeDiff+BMIDiff,study1Over)
summary(fit)
##
## Call:
## lm(formula = AvgAH ~ AgeDiff + BMIDiff, data = study1Over)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.95388 -0.22750 0.09719 0.29277 3.09758
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -15.993961 21.323638 -0.750 0.454
## AgeDiff 0.004988 0.005366 0.930 0.354
## BMIDiff 0.005779 0.016022 0.361 0.719
##
## Residual standard error: 0.8582 on 146 degrees of freedom
## Multiple R-squared: 0.007218, Adjusted R-squared: -0.006382
## F-statistic: 0.5307 on 2 and 146 DF, p-value: 0.5893
vif(fit)
## AgeDiff BMIDiff
## 1.007785 1.007785
shapiro.test(fit$res) ## Shapiro p > 0.01
##
## Shapiro-Wilk normality test
##
## data: fit$res
## W = 0.89753, p-value = 1.049e-08
hist(fit$res)
qqnorm(fit$res)
plot(fit$res)
skewness(fit$res)
## [1] -0.5603742
kurtosis(fit$res)
## [1] 3.148716