####### SETUP #########
library(psych)
library(dplyr)
## Warning: 套件 'dplyr' 是用 R 版本 4.1.3 來建造的
##
## 載入套件：'dplyr'
## 下列物件被遮斷自 'package:stats':
##
## filter, lag
## 下列物件被遮斷自 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyverse)
## Warning: 套件 'tidyverse' 是用 R 版本 4.1.3 來建造的
## -- Attaching packages --------------------------------------- tidyverse 1.3.2 --
## v ggplot2 3.3.6 v purrr 0.3.4
## v tibble 3.1.8 v stringr 1.4.1
## v tidyr 1.2.1 v forcats 0.5.2
## v readr 2.1.2
## Warning: 套件 'ggplot2' 是用 R 版本 4.1.3 來建造的
## Warning: 套件 'tibble' 是用 R 版本 4.1.3 來建造的
## Warning: 套件 'tidyr' 是用 R 版本 4.1.3 來建造的
## Warning: 套件 'readr' 是用 R 版本 4.1.3 來建造的
## Warning: 套件 'purrr' 是用 R 版本 4.1.3 來建造的
## Warning: 套件 'stringr' 是用 R 版本 4.1.3 來建造的
## Warning: 套件 'forcats' 是用 R 版本 4.1.3 來建造的
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x ggplot2::%+%() masks psych::%+%()
## x ggplot2::alpha() masks psych::alpha()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(rstatix)
## Warning: 套件 'rstatix' 是用 R 版本 4.1.3 來建造的
##
## 載入套件：'rstatix'
##
## 下列物件被遮斷自 'package:stats':
##
## filter
library(qualtRics)
## Warning: 套件 'qualtRics' 是用 R 版本 4.1.3 來建造的
library(reshape2)
## Warning: 套件 'reshape2' 是用 R 版本 4.1.3 來建造的
##
## 載入套件：'reshape2'
##
## 下列物件被遮斷自 'package:tidyr':
##
## smiths
library(stringi)
## Warning: 套件 'stringi' 是用 R 版本 4.1.2 來建造的
library(boot)
##
## 載入套件：'boot'
##
## 下列物件被遮斷自 'package:psych':
##
## logit
library(beepr)
## Warning: 套件 'beepr' 是用 R 版本 4.1.3 來建造的
library(multcomp)
## Warning: 套件 'multcomp' 是用 R 版本 4.1.3 來建造的
## 載入需要的套件：mvtnorm
## 載入需要的套件：survival
##
## 載入套件：'survival'
##
## 下列物件被遮斷自 'package:boot':
##
## aml
##
## 載入需要的套件：TH.data
## Warning: 套件 'TH.data' 是用 R 版本 4.1.3 來建造的
## 載入需要的套件：MASS
##
## 載入套件：'MASS'
##
## 下列物件被遮斷自 'package:rstatix':
##
## select
##
## 下列物件被遮斷自 'package:dplyr':
##
## select
##
##
## 載入套件：'TH.data'
##
## 下列物件被遮斷自 'package:MASS':
##
## geyser
library(nlme)
##
## 載入套件：'nlme'
##
## 下列物件被遮斷自 'package:dplyr':
##
## collapse
library(lme4)
## Warning: 套件 'lme4' 是用 R 版本 4.1.3 來建造的
## 載入需要的套件：Matrix
## Warning: 套件 'Matrix' 是用 R 版本 4.1.3 來建造的
##
## 載入套件：'Matrix'
##
## 下列物件被遮斷自 'package:tidyr':
##
## expand, pack, unpack
##
##
## 載入套件：'lme4'
##
## 下列物件被遮斷自 'package:nlme':
##
## lmList
library(emmeans)
## Warning: 套件 'emmeans' 是用 R 版本 4.1.3 來建造的
library(MuMIn)
## Warning: 套件 'MuMIn' 是用 R 版本 4.1.3 來建造的
library("cowplot")
## Warning: 套件 'cowplot' 是用 R 版本 4.1.3 來建造的
library(ggpubr)
## Warning: 套件 'ggpubr' 是用 R 版本 4.1.3 來建造的
##
## 載入套件：'ggpubr'
##
## 下列物件被遮斷自 'package:cowplot':
##
## get_legend
library(factoextra)
## Warning: 套件 'factoextra' 是用 R 版本 4.1.3 來建造的
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(cluster)
library(fossil)
## Warning: 套件 'fossil' 是用 R 版本 4.1.3 來建造的
## 載入需要的套件：sp
## Warning: 套件 'sp' 是用 R 版本 4.1.3 來建造的
## 載入需要的套件：maps
## Warning: 套件 'maps' 是用 R 版本 4.1.3 來建造的
##
## 載入套件：'maps'
##
## 下列物件被遮斷自 'package:cluster':
##
## votes.repub
##
## 下列物件被遮斷自 'package:purrr':
##
## map
##
## 載入需要的套件：shapefiles
## Warning: 套件 'shapefiles' 是用 R 版本 4.1.3 來建造的
## 載入需要的套件：foreign
##
## 載入套件：'shapefiles'
##
## 下列物件被遮斷自 'package:foreign':
##
## read.dbf, write.dbf
##
##
## 載入套件：'fossil'
##
## 下列物件被遮斷自 'package:psych':
##
## manhattan
library(mclust)
## Warning: 套件 'mclust' 是用 R 版本 4.1.3 來建造的
## Package 'mclust' version 6.0.0
## Type 'citation("mclust")' for citing this R package in publications.
##
## 載入套件：'mclust'
##
## 下列物件被遮斷自 'package:maps':
##
## map
##
## 下列物件被遮斷自 'package:mvtnorm':
##
## dmvnorm
##
## 下列物件被遮斷自 'package:purrr':
##
## map
##
## 下列物件被遮斷自 'package:psych':
##
## sim
## Import data ##
##cleanandcomplete.csv
longall <- read.csv("C:/Users/Li jing-yi/Desktop/R/cleanedandcomplete.csv")
longall$country <- as.factor(longall$country)
longall$recipient <- as.factor(longall$recipient)
longall$act <- as.factor(longall$act)
describe(longall)
## vars n mean sd median trimmed mad min max range
## act* 1 4676 292.77 168.78 293.00 292.75 216.46 1.00 585.00 584.00
## country* 2 4680 1.50 0.50 1.50 1.50 0.74 1.00 2.00 1.00
## recipient* 3 4680 2.50 1.12 2.50 2.50 1.48 1.00 4.00 3.00
## kindness 4 3385 6.44 0.95 6.45 6.46 1.08 3.37 8.55 5.17
## cost 5 3383 2.32 1.00 2.23 2.25 0.95 0.31 7.05 6.73
## benefit 6 3385 6.41 1.01 6.46 6.45 1.04 2.08 8.81 6.73
## skew kurtosis se
## act* 0.00 -1.20 2.47
## country* 0.00 -2.00 0.01
## recipient* 0.00 -1.36 0.02
## kindness -0.23 -0.60 0.02
## cost 0.80 1.19 0.02
## benefit -0.43 0.06 0.02
describeBy(longall$benefit,
list(longall$recipient, longall$country))
##
## Descriptive statistics by group
## : colleague
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 327 5.66 0.78 5.76 5.71 0.78 2.11 7.52 5.41 -0.78 1.25 0.04
## ------------------------------------------------------------
## : family
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 536 5.84 0.76 5.93 5.9 0.69 2.72 7.43 4.72 -0.73 0.72 0.03
## ------------------------------------------------------------
## : friend
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 534 5.98 0.78 6.08 6.04 0.77 2.98 7.64 4.65 -0.75 0.61 0.03
## ------------------------------------------------------------
## : stranger
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 295 5.48 0.91 5.53 5.53 0.92 2.08 7.37 5.29 -0.53 0.31 0.05
## ------------------------------------------------------------
## : colleague
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 327 6.78 0.83 6.9 6.83 0.8 3.51 8.45 4.94 -0.72 0.57 0.05
## ------------------------------------------------------------
## : family
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 536 7.12 0.69 7.19 7.17 0.67 4.37 8.61 4.24 -0.75 0.74 0.03
## ------------------------------------------------------------
## : friend
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 534 7.21 0.7 7.3 7.26 0.62 4.19 8.81 4.62 -0.79 1.11 0.03
## ------------------------------------------------------------
## : stranger
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 296 6.82 0.86 6.89 6.87 0.81 3.62 8.56 4.94 -0.65 0.58 0.05
describeBy(longall$benefit,
list(longall$recipient))
##
## Descriptive statistics by group
## : colleague
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 654 6.22 0.98 6.23 6.25 1 2.11 8.45 6.34 -0.34 0.13 0.04
## ------------------------------------------------------------
## : family
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1072 6.48 0.97 6.51 6.52 1.01 2.72 8.61 5.89 -0.41 -0.04 0.03
## ------------------------------------------------------------
## : friend
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1068 6.6 0.96 6.65 6.64 0.99 2.98 8.81 5.83 -0.48 0.09 0.03
## ------------------------------------------------------------
## : stranger
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 591 6.15 1.11 6.26 6.19 1.17 2.08 8.56 6.48 -0.36 -0.1 0.05
describeBy(longall$benefit,
list(longall$country))
##
## Descriptive statistics by group
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1692 5.79 0.82 5.91 5.84 0.8 2.08 7.64 5.56 -0.72 0.74 0.02
## ------------------------------------------------------------
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1693 7.03 0.77 7.14 7.09 0.72 3.51 8.81 5.3 -0.81 0.97 0.02
describeBy(longall$cost,
list(longall$recipient, longall$country))
##
## Descriptive statistics by group
## : colleague
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 327 1.76 0.74 1.63 1.69 0.7 0.62 5.1 4.47 1.17 2.18 0.04
## ------------------------------------------------------------
## : family
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 535 1.63 0.81 1.46 1.53 0.7 0.35 4.88 4.53 1.22 1.63 0.03
## ------------------------------------------------------------
## : friend
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 534 1.71 0.84 1.52 1.61 0.7 0.31 5.02 4.7 1.19 1.53 0.04
## ------------------------------------------------------------
## : stranger
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 295 2.28 0.85 2.17 2.19 0.81 0.82 5.58 4.76 1.19 2.15 0.05
## ------------------------------------------------------------
## : colleague
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 327 2.89 0.77 2.8 2.83 0.74 1.46 7.05 5.59 1.32 4.12 0.04
## ------------------------------------------------------------
## : family
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 535 2.8 0.85 2.63 2.69 0.72 1.54 6.4 4.86 1.25 1.63 0.04
## ------------------------------------------------------------
## : friend
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 534 2.76 0.87 2.55 2.65 0.72 1.41 6.77 5.37 1.5 3.16 0.04
## ------------------------------------------------------------
## : stranger
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 296 3.01 0.94 2.9 2.92 0.8 1.54 6.7 5.16 1.2 2.09 0.05
describeBy(longall$cost,
list(longall$recipient))
##
## Descriptive statistics by group
## : colleague
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 654 2.33 0.94 2.29 2.28 0.97 0.62 7.05 6.42 0.69 1.26 0.04
## ------------------------------------------------------------
## : family
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1070 2.21 1.01 2.13 2.14 0.98 0.35 6.4 6.05 0.75 0.68 0.03
## ------------------------------------------------------------
## : friend
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1068 2.24 1 2.15 2.16 0.95 0.31 6.77 6.46 0.9 1.46 0.03
## ------------------------------------------------------------
## : stranger
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 591 2.65 0.97 2.54 2.56 0.88 0.82 6.7 5.88 1.05 1.78 0.04
describeBy(longall$cost,
list(longall$country))
##
## Descriptive statistics by group
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1691 1.8 0.84 1.63 1.71 0.77 0.31 5.58 5.26 1.13 1.66 0.02
## ------------------------------------------------------------
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1692 2.84 0.87 2.67 2.74 0.74 1.41 7.05 5.64 1.33 2.6 0.02
describeBy(longall$kindness,
list(longall$recipient, longall$country))
##
## Descriptive statistics by group
## : colleague
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 327 5.68 0.6 5.76 5.71 0.56 3.37 7.08 3.71 -0.51 0.17 0.03
## ------------------------------------------------------------
## : family
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 536 5.71 0.61 5.76 5.74 0.59 3.98 7.12 3.14 -0.37 -0.21 0.03
## ------------------------------------------------------------
## : friend
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 534 5.88 0.65 5.93 5.91 0.66 3.94 7.61 3.67 -0.37 -0.1 0.03
## ------------------------------------------------------------
## : stranger
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 295 5.61 0.75 5.69 5.65 0.72 3.41 7.19 3.79 -0.46 -0.08 0.04
## ------------------------------------------------------------
## : colleague
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 327 7 0.59 7.09 7.04 0.54 4.03 8.21 4.18 -0.9 1.8 0.03
## ------------------------------------------------------------
## : family
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 536 7.22 0.58 7.29 7.26 0.57 4.19 8.55 4.35 -0.77 1.39 0.03
## ------------------------------------------------------------
## : friend
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 534 7.25 0.57 7.31 7.29 0.56 4.25 8.42 4.17 -0.8 1.4 0.02
## ------------------------------------------------------------
## : stranger
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 296 6.88 0.77 6.93 6.94 0.76 4.39 8.53 4.14 -0.58 0.1 0.04
describeBy(longall$kindness,
list(longall$recipient))
##
## Descriptive statistics by group
## : colleague
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 654 6.34 0.89 6.36 6.37 1.02 3.37 8.21 4.84 -0.24 -0.57 0.03
## ------------------------------------------------------------
## : family
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1072 6.47 0.96 6.45 6.49 1.17 3.98 8.55 4.56 -0.17 -0.81 0.03
## ------------------------------------------------------------
## : friend
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1068 6.57 0.92 6.6 6.6 1.04 3.94 8.42 4.48 -0.29 -0.62 0.03
## ------------------------------------------------------------
## : stranger
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 591 6.25 0.99 6.31 6.27 1.06 3.41 8.53 5.12 -0.21 -0.42 0.04
describeBy(longall$kindness,
list(longall$country))
##
## Descriptive statistics by group
## : UK
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1692 5.74 0.65 5.79 5.77 0.64 3.37 7.61 4.24 -0.43 0.07 0.02
## ------------------------------------------------------------
## : USA
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1693 7.13 0.63 7.21 7.17 0.59 4.03 8.55 4.51 -0.84 1.3 0.02
overkind <-longall %>%
arrange(desc(kindness)) %>%
slice(1:11)
overkindrec <-longall %>%
arrange(desc(kindness)) %>%
group_by(recipient) %>%
slice(1:11)
##Intercept Only Model##
iomodel <- gls(kindness ~ 1,
data=longall,
method = "ML",
na.action = "na.omit")
summary(iomodel)
## Generalized least squares fit by maximum likelihood
## Model: kindness ~ 1
## Data: longall
## AIC BIC logLik
## 9230.016 9242.27 -4613.008
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 6.43504 0.01625156 395.9643 0
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.23958517 -0.72866614 0.01461647 0.82369710 2.23381599
##
## Residual standard error: 0.9453886
## Degrees of freedom: 3385 total; 3384 residual
##Random intercept only model (recipients)##
riomodel <- lme(kindness ~ 1,
data=longall,
method = "ML",
na.action = "na.omit",
random = ~1|recipient)
summary(riomodel)
## Linear mixed-effects model fit by maximum likelihood
## Data: longall
## AIC BIC logLik
## 9195.051 9213.433 -4594.526
##
## Random effects:
## Formula: ~1 | recipient
## (Intercept) Residual
## StdDev: 0.1155086 0.9388009
##
## Fixed effects: kindness ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 6.407388 0.06012858 3381 106.5614 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.16774937 -0.71959031 0.01775967 0.82842058 2.41084012
##
## Number of Observations: 3385
## Number of Groups: 4
anova(iomodel, riomodel)
## Model df AIC BIC logLik Test L.Ratio p-value
## iomodel 1 2 9230.016 9242.270 -4613.008
## riomodel 2 3 9195.051 9213.433 -4594.526 1 vs 2 36.96442 <.0001
r.squaredGLMM(riomodel)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
## R2m R2c
## [1,] 0 0.01491271
#followups for significant recipient effect
recim <- lmer(kindness~factor(recipient)+(1|act), data=longall)
anova(recim)
## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## factor(recipient) 3 43.68 14.56 21.585
rec.emm <- emmeans(recim, "recipient", pbkrtest.limit = 3392)
rec.emm.df <- as.data.frame(summary(rec.emm))
contrast(rec.emm, 'tukey')
## contrast estimate SE df t.ratio p.value
## colleague - family -0.124 0.0423 3090 -2.936 0.0176
## colleague - friend -0.226 0.0423 3079 -5.339 <.0001
## colleague - stranger 0.101 0.0494 3194 2.040 0.1737
## family - friend -0.102 0.0356 2836 -2.855 0.0225
## family - stranger 0.225 0.0441 3159 5.097 <.0001
## friend - stranger 0.327 0.0442 3156 7.397 <.0001
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates
rec.emm.df %>%
ggplot(aes(recipient, emmean, ymin=lower.CL, ymax=upper.CL)) +
geom_pointrange() +
ylab("kindness")
describe.by(longall, group = "recipient")
## Warning: describe.by is deprecated. Please use the describeBy function
##
## Descriptive statistics by group
## recipient: colleague
## vars n mean sd median trimmed mad min max range
## act* 1 1169 292.77 168.84 293.00 292.75 216.46 1.00 585.00 584.00
## country* 2 1170 1.50 0.50 1.50 1.50 0.74 1.00 2.00 1.00
## recipient* 3 1170 1.00 0.00 1.00 1.00 0.00 1.00 1.00 0.00
## kindness 4 654 6.34 0.89 6.36 6.37 1.02 3.37 8.21 4.84
## cost 5 654 2.33 0.94 2.29 2.28 0.97 0.62 7.05 6.42
## benefit 6 654 6.22 0.98 6.23 6.25 1.00 2.11 8.45 6.34
## skew kurtosis se
## act* 0.00 -1.20 4.94
## country* 0.00 -2.00 0.01
## recipient* NaN NaN 0.00
## kindness -0.24 -0.57 0.03
## cost 0.69 1.26 0.04
## benefit -0.34 0.13 0.04
## ------------------------------------------------------------
## recipient: family
## vars n mean sd median trimmed mad min max range
## act* 1 1169 292.77 168.84 293.00 292.75 216.46 1.00 585.00 584.00
## country* 2 1170 1.50 0.50 1.50 1.50 0.74 1.00 2.00 1.00
## recipient* 3 1170 2.00 0.00 2.00 2.00 0.00 2.00 2.00 0.00
## kindness 4 1072 6.47 0.96 6.45 6.49 1.17 3.98 8.55 4.56
## cost 5 1070 2.21 1.01 2.13 2.14 0.98 0.35 6.40 6.05
## benefit 6 1072 6.48 0.97 6.51 6.52 1.01 2.72 8.61 5.89
## skew kurtosis se
## act* 0.00 -1.20 4.94
## country* 0.00 -2.00 0.01
## recipient* NaN NaN 0.00
## kindness -0.17 -0.81 0.03
## cost 0.75 0.68 0.03
## benefit -0.41 -0.04 0.03
## ------------------------------------------------------------
## recipient: friend
## vars n mean sd median trimmed mad min max range
## act* 1 1169 292.77 168.84 293.00 292.75 216.46 1.00 585.00 584.00
## country* 2 1170 1.50 0.50 1.50 1.50 0.74 1.00 2.00 1.00
## recipient* 3 1170 3.00 0.00 3.00 3.00 0.00 3.00 3.00 0.00
## kindness 4 1068 6.57 0.92 6.60 6.60 1.04 3.94 8.42 4.48
## cost 5 1068 2.24 1.00 2.15 2.16 0.95 0.31 6.77 6.46
## benefit 6 1068 6.60 0.96 6.65 6.64 0.99 2.98 8.81 5.83
## skew kurtosis se
## act* 0.00 -1.20 4.94
## country* 0.00 -2.00 0.01
## recipient* NaN NaN 0.00
## kindness -0.29 -0.62 0.03
## cost 0.90 1.46 0.03
## benefit -0.48 0.09 0.03
## ------------------------------------------------------------
## recipient: stranger
## vars n mean sd median trimmed mad min max range
## act* 1 1169 292.77 168.84 293.00 292.75 216.46 1.00 585.00 584.00
## country* 2 1170 1.50 0.50 1.50 1.50 0.74 1.00 2.00 1.00
## recipient* 3 1170 4.00 0.00 4.00 4.00 0.00 4.00 4.00 0.00
## kindness 4 591 6.25 0.99 6.31 6.27 1.06 3.41 8.53 5.12
## cost 5 591 2.65 0.97 2.54 2.56 0.88 0.82 6.70 5.88
## benefit 6 591 6.15 1.11 6.26 6.19 1.17 2.08 8.56 6.48
## skew kurtosis se
## act* 0.00 -1.20 4.94
## country* 0.00 -2.00 0.01
## recipient* NaN NaN 0.00
## kindness -0.21 -0.42 0.04
## cost 1.05 1.78 0.04
## benefit -0.36 -0.10 0.05
##Random intercept only model (country+ recipients)##
riomodel2 <- lme(kindness ~ 1,
data=longall,
method = "ML",
na.action = "na.omit",
random = list (~1|country, ~1|recipient))
summary(riomodel2)
## Linear mixed-effects model fit by maximum likelihood
## Data: longall
## AIC BIC logLik
## 6532.232 6556.74 -3262.116
##
## Random effects:
## Formula: ~1 | country
## (Intercept)
## StdDev: 0.6796816
##
## Formula: ~1 | recipient %in% country
## (Intercept) Residual
## StdDev: 0.1455203 0.6311312
##
## Fixed effects: kindness ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 6.406145 0.4835555 3377 13.248 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.7895889 -0.5962011 0.0982295 0.7027948 2.7577684
##
## Number of Observations: 3385
## Number of Groups:
## country recipient %in% country
## 2 8
anova(riomodel, riomodel2)
## Model df AIC BIC logLik Test L.Ratio p-value
## riomodel 1 3 9195.051 9213.433 -4594.526
## riomodel2 2 4 6532.232 6556.740 -3262.116 1 vs 2 2664.819 <.0001
r.squaredGLMM(riomodel2)
## R2m R2c
## [1,] 0 0.5481109
#followups for significant country effect
longall$country <-as.numeric(longall$country)
counm <- lmer(kindness~factor(country)+ (1|recipient), data=longall)
anova(counm)
## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## factor(country) 1 1627.6 1627.6 4068.5
coun.emm <- emmeans(counm, "country",pbkrtest.limit = 3392)
coun.emm.df <- as.data.frame(summary(coun.emm))
coun.emm.df$country <- as.factor(coun.emm.df$country)
coun.emm.df %>%
ggplot(aes(country, emmean, ymin=lower.CL, ymax=upper.CL)) +
geom_pointrange() +
ylab("kindness")
longall$country <-as.factor(longall$country)
longall$country <- dplyr::recode(longall$country, "1" = "USA",
"2" = "UK")
describe.by(longall, group = "country")
## Warning: describe.by is deprecated. Please use the describeBy function
##
## Descriptive statistics by group
## country: USA
## vars n mean sd median trimmed mad min max range
## act* 1 2336 292.54 168.69 292.50 292.50 216.46 1.00 585.00 584.00
## country* 2 2340 1.00 0.00 1.00 1.00 0.00 1.00 1.00 0.00
## recipient* 3 2340 2.50 1.12 2.50 2.50 1.48 1.00 4.00 3.00
## kindness 4 1692 5.74 0.65 5.79 5.77 0.64 3.37 7.61 4.24
## cost 5 1691 1.80 0.84 1.63 1.71 0.77 0.31 5.58 5.26
## benefit 6 1692 5.79 0.82 5.91 5.84 0.80 2.08 7.64 5.56
## skew kurtosis se
## act* 0.00 -1.20 3.49
## country* NaN NaN 0.00
## recipient* 0.00 -1.36 0.02
## kindness -0.43 0.07 0.02
## cost 1.13 1.66 0.02
## benefit -0.72 0.74 0.02
## ------------------------------------------------------------
## country: UK
## vars n mean sd median trimmed mad min max range
## act* 1 2340 293.00 168.91 293.00 293.00 216.46 1.00 585.00 584.00
## country* 2 2340 2.00 0.00 2.00 2.00 0.00 2.00 2.00 0.00
## recipient* 3 2340 2.50 1.12 2.50 2.50 1.48 1.00 4.00 3.00
## kindness 4 1693 7.13 0.63 7.21 7.17 0.59 4.03 8.55 4.51
## cost 5 1692 2.84 0.87 2.67 2.74 0.74 1.41 7.05 5.64
## benefit 6 1693 7.03 0.77 7.14 7.09 0.72 3.51 8.81 5.30
## skew kurtosis se
## act* 0.00 -1.20 3.49
## country* NaN NaN 0.00
## recipient* 0.00 -1.36 0.02
## kindness -0.84 1.30 0.02
## cost 1.33 2.60 0.02
## benefit -0.81 0.97 0.02
####Random intercept only model (acts+ recipients + country)####
riomodel3 <- lme(kindness ~ 1,
data=longall,
method = "ML",
na.action = "na.omit",
random = list (~1|act, ~1|country, ~1|recipient))
summary(riomodel3)
## Linear mixed-effects model fit by maximum likelihood
## Data: longall
## AIC BIC logLik
## 5592.325 5622.961 -2791.163
##
## Random effects:
## Formula: ~1 | act
## (Intercept)
## StdDev: 0.0001355577
##
## Formula: ~1 | country %in% act
## (Intercept)
## StdDev: 0.8927543
##
## Formula: ~1 | recipient %in% country %in% act
## (Intercept) Residual
## StdDev: 0.3173234 0.06788708
##
## Fixed effects: kindness ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 6.415831 0.02678934 2216 239.492 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.140210288 -0.099430958 0.009527817 0.110826598 0.572600573
##
## Number of Observations: 3385
## Number of Groups:
## act country %in% act
## 585 1169
## recipient %in% country %in% act
## 3385
anova(riomodel3, riomodel2)
## Model df AIC BIC logLik Test L.Ratio p-value
## riomodel3 1 5 5592.325 5622.961 -2791.163
## riomodel2 2 4 6532.232 6556.740 -3262.116 1 vs 2 941.9065 <.0001
r.squaredGLMM(riomodel3)
## R2m R2c
## [1,] 0 0.9948924
##Predictor model (benefit)##
premodel1 <- lme(kindness ~ benefit,
data=longall,
method = "ML",
na.action = "na.omit",
random = list (~1|act, ~1|country, ~1|recipient))
summary(premodel1)
## Linear mixed-effects model fit by maximum likelihood
## Data: longall
## AIC BIC logLik
## 4001.953 4038.715 -1994.976
##
## Random effects:
## Formula: ~1 | act
## (Intercept)
## StdDev: 9.83218e-05
##
## Formula: ~1 | country %in% act
## (Intercept)
## StdDev: 0.5134714
##
## Formula: ~1 | recipient %in% country %in% act
## (Intercept) Residual
## StdDev: 0.2968044 0.02049741
##
## Fixed effects: kindness ~ benefit
## Value Std.Error DF t-value p-value
## (Intercept) 3.198699 0.06593517 2215 48.51279 0
## benefit 0.505523 0.01004399 2215 50.33084 0
## Correlation:
## (Intr)
## benefit -0.97
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.246315897 -0.034485722 0.002494742 0.038633540 0.215291560
##
## Number of Observations: 3385
## Number of Groups:
## act country %in% act
## 585 1169
## recipient %in% country %in% act
## 3385
anova(riomodel3, premodel1)
## Model df AIC BIC logLik Test L.Ratio p-value
## riomodel3 1 5 5592.325 5622.961 -2791.163
## premodel1 2 6 4001.953 4038.715 -1994.976 1 vs 2 1592.373 <.0001
plot.lme(premodel1, kindness ~ benefit | recipient + country)
r.squaredGLMM(premodel1)
## R2m R2c
## [1,] 0.425002 0.999314
library(Matrix)
library(lme4)
##Predictor model (benefit)##
premodel1 <- lme(kindness ~ benefit,
data=longall,
method = "ML",
na.action = "na.omit",
random = list (~1|act, ~1|country, ~1|recipient))
summary(premodel1)
## Linear mixed-effects model fit by maximum likelihood
## Data: longall
## AIC BIC logLik
## 4001.953 4038.715 -1994.976
##
## Random effects:
## Formula: ~1 | act
## (Intercept)
## StdDev: 9.83218e-05
##
## Formula: ~1 | country %in% act
## (Intercept)
## StdDev: 0.5134714
##
## Formula: ~1 | recipient %in% country %in% act
## (Intercept) Residual
## StdDev: 0.2968044 0.02049741
##
## Fixed effects: kindness ~ benefit
## Value Std.Error DF t-value p-value
## (Intercept) 3.198699 0.06593517 2215 48.51279 0
## benefit 0.505523 0.01004399 2215 50.33084 0
## Correlation:
## (Intr)
## benefit -0.97
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.246315897 -0.034485722 0.002494742 0.038633540 0.215291560
##
## Number of Observations: 3385
## Number of Groups:
## act country %in% act
## 585 1169
## recipient %in% country %in% act
## 3385
anova(riomodel3, premodel1)
## Model df AIC BIC logLik Test L.Ratio p-value
## riomodel3 1 5 5592.325 5622.961 -2791.163
## premodel1 2 6 4001.953 4038.715 -1994.976 1 vs 2 1592.373 <.0001
plot.lme(premodel1, kindness ~ benefit | recipient + country)
r.squaredGLMM(premodel1)
## R2m R2c
## [1,] 0.425002 0.999314
##Predictor model (benefit + cost)##
premodel2 <- lme(kindness ~ benefit + cost,
data=longall,
method = "ML",
na.action = "na.omit",
random = list (~1|act, ~1|country, ~1|recipient))
summary(premodel2)
## Linear mixed-effects model fit by maximum likelihood
## Data: longall
## AIC BIC logLik
## 3605.84 3648.726 -1795.92
##
## Random effects:
## Formula: ~1 | act
## (Intercept)
## StdDev: 0.2227582
##
## Formula: ~1 | country %in% act
## (Intercept)
## StdDev: 0.22953
##
## Formula: ~1 | recipient %in% country %in% act
## (Intercept) Residual
## StdDev: 0.3084346 0.129045
##
## Fixed effects: kindness ~ benefit + cost
## Value Std.Error DF t-value p-value
## (Intercept) 1.6933018 0.05862660 2212 28.88283 0
## benefit 0.6391430 0.00892854 2212 71.58428 0
## cost 0.2772236 0.00993388 2212 27.90688 0
## Correlation:
## (Intr) beneft
## benefit -0.893
## cost -0.211 -0.194
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.72015001 -0.19813535 0.01264534 0.23018399 1.57964695
##
## Number of Observations: 3383
## Number of Groups:
## act country %in% act
## 585 1169
## recipient %in% country %in% act
## 3383
plot.lme(premodel2, kindness ~ cost | recipient + country)
r.squaredGLMM(premodel2)
## R2m R2c
## [1,] 0.7323825 0.9791838
########Predictor model (benefit + cost + interaction)######
##Variable centering##
longall$cbenefit <- longall$benefit - mean(longall$benefit, na.rm=TRUE)
longall$ccost <- longall$cost - mean(longall$cost, na.rm=TRUE)
##Model##
premodel3 <- lme(kindness ~ cbenefit * ccost,
data=longall,
method = "ML",
na.action = "na.omit",
random = list (~1|act, ~1|country, ~1|recipient))
summary(premodel3)
## Linear mixed-effects model fit by maximum likelihood
## Data: longall
## AIC BIC logLik
## 3597.932 3646.944 -1790.966
##
## Random effects:
## Formula: ~1 | act
## (Intercept)
## StdDev: 0.2108745
##
## Formula: ~1 | country %in% act
## (Intercept)
## StdDev: 0.2536281
##
## Formula: ~1 | recipient %in% country %in% act
## (Intercept) Residual
## StdDev: 0.3025844 0.1292639
##
## Fixed effects: kindness ~ cbenefit * ccost
## Value Std.Error DF t-value p-value
## (Intercept) 6.424254 0.013127044 2211 489.3908 0e+00
## cbenefit 0.625788 0.009043996 2211 69.1937 0e+00
## ccost 0.267970 0.009999618 2211 26.7981 0e+00
## cbenefit:ccost 0.030050 0.008820540 2211 3.4068 7e-04
## Correlation:
## (Intr) cbenft ccost
## cbenefit 0.046
## ccost -0.035 -0.159
## cbenefit:ccost -0.175 -0.102 0.004
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.66543107 -0.20426369 0.00867212 0.23214056 1.48885044
##
## Number of Observations: 3383
## Number of Groups:
## act country %in% act
## 585 1169
## recipient %in% country %in% act
## 3383
anova(premodel3, premodel2)
## Model df AIC BIC logLik Test L.Ratio p-value
## premodel3 1 8 3597.932 3646.944 -1790.966
## premodel2 2 7 3605.840 3648.726 -1795.920 1 vs 2 9.908651 0.0016
plot.lme(premodel3, kindness ~ benefit*cost | recipient + country)
r.squaredGLMM(premodel3)
## R2m R2c
## [1,] 0.7227575 0.9786582
rsmodel1a <- lme(kindness ~ cbenefit*ccost,
data=longall,
method = "ML",
na.action = "na.omit",
random = list (~1|country, ~1|recipient),
control =list(msMaxIter = 1000, msMaxEval = 1000))
summary(rsmodel1a)
## Linear mixed-effects model fit by maximum likelihood
## Data: longall
## AIC BIC logLik
## 3809.475 3852.361 -1897.737
##
## Random effects:
## Formula: ~1 | country
## (Intercept)
## StdDev: 0.2261358
##
## Formula: ~1 | recipient %in% country
## (Intercept) Residual
## StdDev: 0.06778803 0.4223298
##
## Fixed effects: kindness ~ cbenefit * ccost
## Value Std.Error DF t-value p-value
## (Intercept) 6.424681 0.16196983 3372 39.66591 0.0000
## cbenefit 0.581693 0.00945142 3372 61.54553 0.0000
## ccost 0.184826 0.00881290 3372 20.97225 0.0000
## cbenefit:ccost 0.008958 0.00782787 3372 1.14435 0.2526
## Correlation:
## (Intr) cbenft ccost
## cbenefit 0.004
## ccost 0.000 0.084
## cbenefit:ccost -0.012 -0.135 -0.162
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.99879995 -0.56177971 0.09229664 0.67242841 4.30055840
##
## Number of Observations: 3383
## Number of Groups:
## country recipient %in% country
## 2 8
rsmodel1 <- lme(kindness ~ cbenefit*ccost,
data=longall,
method = "ML",
na.action = "na.omit",
random = list (~1|country, ~cbenefit|recipient),
control =list(msMaxIter = 1000, msMaxEval = 1000))
summary(rsmodel1)
## Linear mixed-effects model fit by maximum likelihood
## Data: longall
## AIC BIC logLik
## 3701.59 3756.729 -1841.795
##
## Random effects:
## Formula: ~1 | country
## (Intercept)
## StdDev: 0.2082933
##
## Formula: ~cbenefit | recipient %in% country
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.1064273 (Intr)
## cbenefit 0.1113263 0.125
## Residual 0.4136363
##
## Fixed effects: kindness ~ cbenefit * ccost
## Value Std.Error DF t-value p-value
## (Intercept) 6.407340 0.15239861 3372 42.04330 0.0000
## cbenefit 0.591065 0.04050747 3372 14.59152 0.0000
## ccost 0.187229 0.00873502 3372 21.43424 0.0000
## cbenefit:ccost -0.018764 0.00863794 3372 -2.17225 0.0299
## Correlation:
## (Intr) cbenft ccost
## cbenefit 0.031
## ccost -0.001 0.024
## cbenefit:ccost 0.001 -0.046 -0.185
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.25270850 -0.57507736 0.08540641 0.68831853 3.07041273
##
## Number of Observations: 3383
## Number of Groups:
## country recipient %in% country
## 2 8
coef(rsmodel1)
## (Intercept) cbenefit ccost cbenefit:ccost
## USA/colleague 6.177453 0.5139913 0.1872286 -0.0187638
## USA/family 6.065449 0.3788506 0.1872286 -0.0187638
## USA/friend 6.253997 0.6010328 0.1872286 -0.0187638
## USA/stranger 6.252297 0.6841688 0.1872286 -0.0187638
## UK/colleague 6.686442 0.5393286 0.1872286 -0.0187638
## UK/family 6.698749 0.6196598 0.1872286 -0.0187638
## UK/friend 6.667450 0.6291690 0.1872286 -0.0187638
## UK/stranger 6.456885 0.7623214 0.1872286 -0.0187638
rsmodel2 <- lme(kindness ~ cbenefit*ccost,
data=longall,
method = "ML",
na.action = "na.omit",
random = list (~1|country, ~cbenefit + ccost|recipient),
control =list(msMaxIter = 1000, msMaxEval = 1000)
)
summary(rsmodel2)
## Linear mixed-effects model fit by maximum likelihood
## Data: longall
## AIC BIC logLik
## 3641.754 3715.272 -1808.877
##
## Random effects:
## Formula: ~1 | country
## (Intercept)
## StdDev: 0.2404818
##
## Formula: ~cbenefit + ccost | recipient %in% country
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.07556837 (Intr) cbenft
## cbenefit 0.09664184 0.010
## ccost 0.08734430 0.328 -0.739
## Residual 0.40903236
##
## Fixed effects: kindness ~ cbenefit * ccost
## Value Std.Error DF t-value p-value
## (Intercept) 6.444895 0.17254377 3372 37.35223 0.0000
## cbenefit 0.589688 0.03546867 3372 16.62560 0.0000
## ccost 0.185749 0.03215539 3372 5.77660 0.0000
## cbenefit:ccost 0.020233 0.01011199 3372 2.00092 0.0455
## Correlation:
## (Intr) cbenft ccost
## cbenefit 0.002
## ccost 0.048 -0.676
## cbenefit:ccost 0.012 -0.051 -0.037
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.6408421 -0.5629546 0.0853122 0.6852110 3.0782168
##
## Number of Observations: 3383
## Number of Groups:
## country recipient %in% country
## 2 8
anova(rsmodel1, rsmodel2)
## Model df AIC BIC logLik Test L.Ratio p-value
## rsmodel1 1 9 3701.590 3756.729 -1841.795
## rsmodel2 2 12 3641.754 3715.272 -1808.877 1 vs 2 65.83664 <.0001