library(reshape);library(ggplot2);library(dplyr)
##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:reshape':
##
## rename
##
## The following objects are masked from 'package:stats':
##
## filter, lag
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(corrplot);library(Hmisc);library(candisc);library(lm.beta)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
##
## The following objects are masked from 'package:dplyr':
##
## combine, src, summarize
##
## The following objects are masked from 'package:base':
##
## format.pval, round.POSIXt, trunc.POSIXt, units
##
## Loading required package: car
## Loading required package: heplots
##
## Attaching package: 'candisc'
##
## The following object is masked from 'package:stats':
##
## cancor
dta <- read.csv("data/socialdevelop.csv",h=T)
head(dta);tail(dta)
dta <- dta[1:100,1:9]
dta <- na.omit(dta)
dim(dta)
dta$age <- dta$age/365
aggregate(dta[,5:9], by=list(dta$sex), FUN=mean)
## Group.1 ToM_Neutral ToM_Emapthy express_thoughts prosocial_behavior
## 1 1 4.425000 4.700000 7.275000 2.100000
## 2 2 4.439024 4.585366 8.707317 1.634146
## offensive_behavior
## 1 2.050000
## 2 1.268293
aggregate(dta[,5:9], by=list(dta$sex), FUN=sd)
## Group.1 ToM_Neutral ToM_Emapthy express_thoughts prosocial_behavior
## 1 1 2.086434 2.462019 7.299412 2.436896
## 2 2 1.962763 2.179170 6.929083 1.639453
## offensive_behavior
## 1 2.286527
## 2 1.923855
mean(dta$age);sd(dta$age)
## [1] 4.688889
## [1] 0.9112838
dtalong <- melt(dta,id=c(colnames(dta[,1:6])))
head(dtalong)
## id sex age class ToM_Neutral ToM_Emapthy variable value
## 1 1 1 6.268493 big 5 9 express_thoughts 6
## 2 2 2 6.136986 big 7 4 express_thoughts 4
## 3 3 2 6.095890 big 6 0 express_thoughts 7
## 4 5 1 6.084932 big 5 3 express_thoughts 14
## 5 6 1 5.958904 big 9 10 express_thoughts 6
## 6 7 2 5.898630 big 8 7 express_thoughts 10
dtalong <- rename(dtalong,factors = variable, scores=value)
ggplot(dtalong,aes(x=age,y=scores,col=as.factor(sex)))+
geom_point()+
facet_grid(.~factors)+
geom_smooth(method="glm")+
scale_colour_brewer(palette="Set3")
summary(glm(express_thoughts~age+sex,data=dta,family = "poisson"))
##
## Call:
## glm(formula = express_thoughts ~ age + sex, family = "poisson",
## data = dta)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8561 -1.9999 -0.2574 0.8977 6.6075
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.32907 0.24814 1.326 0.185
## age 0.31802 0.04352 7.308 2.71e-13 ***
## sex 0.14164 0.07901 1.793 0.073 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 517.61 on 80 degrees of freedom
## Residual deviance: 458.93 on 78 degrees of freedom
## AIC: 730.42
##
## Number of Fisher Scoring iterations: 5
summary(glm(prosocial_behavior~age+sex,data=dta,family = "poisson"))
##
## Call:
## glm(formula = prosocial_behavior ~ age + sex, family = "poisson",
## data = dta)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1189 -0.9279 -0.4762 0.3159 4.7616
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.80512 0.46786 1.721 0.0853 .
## age 0.04227 0.08849 0.478 0.6328
## sex -0.25833 0.16448 -1.571 0.1163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 153.98 on 80 degrees of freedom
## Residual deviance: 151.40 on 78 degrees of freedom
## AIC: 320.27
##
## Number of Fisher Scoring iterations: 5
summary(glm(offensive_behavior~age+sex,data=dta,family = "poisson"))
##
## Call:
## glm(formula = offensive_behavior ~ age + sex, family = "poisson",
## data = dta)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5801 -1.5785 -0.6494 0.5328 3.4539
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.38831 0.52125 -0.745 0.456292
## age 0.34064 0.09193 3.706 0.000211 ***
## sex -0.51930 0.17731 -2.929 0.003402 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 212.77 on 80 degrees of freedom
## Residual deviance: 191.46 on 78 degrees of freedom
## AIC: 324.31
##
## Number of Fisher Scoring iterations: 5
dtatom <- melt(dta[1:6],id=c(colnames(dta[,1:4])))
aggregate(dta[,5:6], by=list(dta$sex), FUN=mean)
## Group.1 ToM_Neutral ToM_Emapthy
## 1 1 4.425000 4.700000
## 2 2 4.439024 4.585366
aggregate(dta[,5:6], by=list(dta$sex), FUN=sd)
## Group.1 ToM_Neutral ToM_Emapthy
## 1 1 2.086434 2.462019
## 2 2 1.962763 2.179170
head(dtatom)
## id sex age class variable value
## 1 1 1 6.268493 big ToM_Neutral 5
## 2 2 2 6.136986 big ToM_Neutral 7
## 3 3 2 6.095890 big ToM_Neutral 6
## 4 5 1 6.084932 big ToM_Neutral 5
## 5 6 1 5.958904 big ToM_Neutral 9
## 6 7 2 5.898630 big ToM_Neutral 8
dtatom <- rename(dtatom,condition = variable, scores=value)
ggplot(dtatom,aes(x=age,y=scores,col=as.factor(sex)))+
geom_point()+
facet_grid(.~condition)+
geom_smooth(method="lm")+
scale_colour_brewer(palette="Set3")
summary(lm.beta(lm(ToM_Neutral~sex+age,data=dta)))
##
## Call:
## lm(formula = ToM_Neutral ~ sex + age, data = dta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5006 -1.2588 -0.2542 1.2445 4.0498
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -0.89299 0.00000 1.10984 -0.805 0.423
## sex -0.20995 -0.05249 0.38259 -0.549 0.585
## age 1.20312 0.54490 0.21121 5.696 2.07e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.712 on 78 degrees of freedom
## Multiple R-squared: 0.2938, Adjusted R-squared: 0.2757
## F-statistic: 16.22 on 2 and 78 DF, p-value: 1.283e-06
summary(lm.beta(lm(ToM_Emapthy~sex+age,data=dta)))
##
## Call:
## lm(formula = ToM_Emapthy ~ sex + age, data = dta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.9296 -1.2465 -0.2101 0.7034 5.6361
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 0.30822 0.00000 1.38785 0.222 0.824829
## sex -0.30493 -0.06643 0.47843 -0.637 0.525754
## age 1.02221 0.40339 0.26412 3.870 0.000224 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.141 on 78 degrees of freedom
## Multiple R-squared: 0.1616, Adjusted R-squared: 0.1401
## F-statistic: 7.519 on 2 and 78 DF, p-value: 0.001033
ggplot(dtatom,aes(x=age,y=scores,col=condition))+
geom_point()+
#facet_grid(.~)+
geom_smooth(method="lm")+
scale_colour_brewer(palette="Set3")
anova(lm(ToM_Neutral~sex+age,data=dta),lm(ToM_Neutral~sex*age,data=dta))
## Analysis of Variance Table
##
## Model 1: ToM_Neutral ~ sex + age
## Model 2: ToM_Neutral ~ sex * age
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 78 228.72
## 2 77 225.91 1 2.812 0.9585 0.3306
anova(lm(ToM_Emapthy~sex+age,data=dta),lm(ToM_Emapthy~sex*age,data=dta))
## Analysis of Variance Table
##
## Model 1: ToM_Emapthy ~ sex + age
## Model 2: ToM_Emapthy ~ sex * age
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 78 357.67
## 2 77 344.04 1 13.629 3.0505 0.0847 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
scaledta <- dta;scaledta[,5:9] <- scale(scaledta[,5:9])
ToM <- scaledta[,5:6];Social_Behavior <- scaledta[,7:9]
canr1 <- cancor(ToM,Social_Behavior)
canr1
##
## Canonical correlation analysis of:
## 2 X variables: ToM_Neutral, ToM_Emapthy
## with 3 Y variables: express_thoughts, prosocial_behavior, offensive_behavior
##
## CanR CanRSQ Eigen percent cum scree
## 1 0.2903 0.08427 0.09202 69.67 69.67 ******************************
## 2 0.1962 0.03851 0.04006 30.33 100.00 *************
##
## Test of H0: The canonical correlations in the
## current row and all that follow are zero
##
## CanR WilksL F df1 df2 p.value
## 1 0.29029 0.88047 1.6649 6 152 0.13328
## 2 0.19625 0.96149 2
coef(canr1,type="both",standardize = T)
## [[1]]
## Xcan1 Xcan2
## ToM_Neutral 0.4870658 1.01070123
## ToM_Emapthy -1.1216812 -0.02411221
##
## [[2]]
## Ycan1 Ycan2
## express_thoughts -0.02822108 0.8489070
## prosocial_behavior -1.01064787 -0.1785497
## offensive_behavior 0.30626474 0.3970981
scaledta$ToM_Neutral_Re <- lm(ToM_Neutral~age+sex,data=scaledta)$residuals
scaledta$ToM_Emapthy_Re <- lm(ToM_Emapthy~age+sex,data=scaledta)$residuals
scaledta$express_thoughts_Re <- glm(express_thoughts~age+sex,data=scaledta)$residuals
scaledta$prosocial_behavior_Re <- lm(prosocial_behavior~age+sex,data=scaledta)$residuals
scaledta$offensive_behavior_Re <- lm(offensive_behavior~age+sex,data=scaledta)$residuals
ToM_re <- scaledta[,10:11];Social_Behavior_re <- scaledta[,12:14]
canr2 <- cancor(ToM_re,Social_Behavior_re)
canr2
##
## Canonical correlation analysis of:
## 2 X variables: ToM_Neutral_Re, ToM_Emapthy_Re
## with 3 Y variables: express_thoughts_Re, prosocial_behavior_Re, offensive_behavior_Re
##
## CanR CanRSQ Eigen percent cum
## 1 0.31772 0.1009489 0.1122839 99.7364 99.74
## 2 0.01722 0.0002966 0.0002967 0.2636 100.00
## scree
## 1 ******************************
## 2
##
## Test of H0: The canonical correlations in the
## current row and all that follow are zero
##
## CanR WilksL F df1 df2 p.value
## 1 0.31772 0.89878 1.3884 6 152 0.22271
## 2 0.01722 0.99970 2
coef(canr2,type="both",standardize = T)
## [[1]]
## Xcan1 Xcan2
## ToM_Neutral_Re 0.2206894 1.0271939
## ToM_Emapthy_Re -1.0453751 -0.1049856
##
## [[2]]
## Ycan1 Ycan2
## express_thoughts_Re 0.1656196 0.666962
## prosocial_behavior_Re -1.0117864 -0.327679
## offensive_behavior_Re 0.4262238 -0.904028