# Alana Survey Data
library(Hmisc)
library(foreign)
library(faux)
library(lme4)
library(nlme)
library(tableone)
library(Epi)
# Set working directory
## setwd('C:/Users/25105517/OneDrive - NSW Health Department/Data/')
## setwd('~/OneDrive - NSW Health Department/Data/')
setwd('~/OneDrive - University of Wollongong/Data/')
## setwd("C:/Users/25105517/OneDrive - University of Wollongong/Data/")
# Load CSV file from working directory
dta <- read.csv('alanaFSurverySept2023.csv', header = T, sep = ',')
# Calculate age from date of survey and DOB
dta$age <- (dta$formdate - dta$DOB)/(365.25*24*60*60)
min(dta$age)## [1] 5.045859max(dta$age)## [1] 12.11773# Create id
dta$id <- 1:length(dta[,1])
# Subset data for analysis
names(dta)## [1] "X" "Q1" "Q10" "Q14" "Q15"
## [6] "Q19" "Q1A" "Q1B" "Q1C" "Q2"
## [11] "Q2_3_TEXT" "Q2_3_0" "Q2_3_1" "Q2_3_2" "Q2_3_3"
## [16] "Q2_3_4" "Q2_3_5" "Q2_3_6" "Q3" "Q30"
## [21] "Q31" "Q32" "Q33" "Q34" "Q35"
## [26] "Q36" "Q4" "Q40" "Q41" "Q42"
## [31] "Q43" "Q44" "Q45" "Q46" "Q5"
## [36] "Q50" "Q51" "Q52" "Q53" "Q54"
## [41] "Q55" "Q56" "Q6" "Q60" "Q61"
## [46] "Q62" "Q63" "Q64" "Q65" "Q66"
## [51] "Q7" "Q8" "Q8_3_TEXT" "Q8_3_0" "Q8_3_1"
## [56] "Q8_3_2" "Q8_3_3" "Q8_3_4" "Q8_3_5" "Q8_3_6"
## [61] "DOB" "Q11" "Q12" "Q12_1_TEXT" "Q12_10"
## [66] "Q12_11" "Q12_12" "Q12_13" "Q12_14" "Q12_15"
## [71] "Q12_16" "Q13" "Q130" "Q131" "Q132"
## [76] "Q133" "Q134" "Q135" "Q136" "Q14_1"
## [81] "Q14_2" "Q14_3" "Q14_4" "Q14_5" "Q14_6"
## [86] "Q14_7" "Q14_8" "Q14_8_TEXT" "Q14_80" "Q14_81"
## [91] "Q14_82" "Q14_83" "Q14_84" "Q14_85" "Q14_86"
## [96] "Q15_1" "Q15_2" "Q15_3" "Q15_4" "Q15_5"
## [101] "Q15_6" "Q16" "Q160" "Q161" "Q162"
## [106] "Q163" "Q164" "Q165" "Q166" "Q17"
## [111] "Q18" "Q180" "Q181" "Q182" "Q183"
## [116] "Q184" "Q185" "Q186" "Q10_1" "Q10_2"
## [121] "Q10_3" "Q10_4" "Q10_5" "Q10_6" "Q10_6_TEXT"
## [126] "Q10_60" "Q10_61" "Q10_62" "Q10_63" "Q10_64"
## [131] "Q10_65" "Q10_66" "Group" "Child.Tray.1" "Child.Tray.2"
## [136] "Child.Tray.3" "Adult.Tray.1" "Adult.Tray.2" "Adult.Tray.3" "formdate"
## [141] "age" "id"describe(dta[,c(141,133,134:139)])## dta[, c(141, 133, 134:139)]
##
## 8 Variables 49 Observations
## --------------------------------------------------------------------------------
## age
## n missing distinct Info Mean Gmd .05 .10
## 49 0 48 1 8.42 2.276 5.288 5.560
## .25 .50 .75 .90 .95
## 6.658 8.739 10.034 10.737 11.114
##
## lowest : 5.045859 5.221081 5.259411 5.330595 5.492129
## highest: 10.951403 11.025325 11.173169 11.603012 12.117728
## --------------------------------------------------------------------------------
## Group
## n missing distinct Info Mean Gmd
## 49 0 3 0.889 1.98 0.9116
##
## Value 1 2 3
## Frequency 17 16 16
## Proportion 0.347 0.327 0.327
## --------------------------------------------------------------------------------
## Child.Tray.1
## n missing distinct
## 49 0 4
##
## Value cookies fruit muesli bar yoghurt
## Frequency 21 10 9 9
## Proportion 0.429 0.204 0.184 0.184
## --------------------------------------------------------------------------------
## Child.Tray.2
## n missing distinct
## 49 0 5
##
## lowest : cookies fruit muesli bar yoghurt
## highest: cookies fruit muesli bar yoghurt
##
## Value cookies fruit muesli bar yoghurt
## Frequency 17 8 16 7 1
## Proportion 0.347 0.163 0.327 0.143 0.020
## --------------------------------------------------------------------------------
## Child.Tray.3
## n missing distinct
## 49 0 4
##
## Value cookies fruit muesli bar
## Frequency 17 22 3 7
## Proportion 0.347 0.449 0.061 0.143
## --------------------------------------------------------------------------------
## Adult.Tray.1
## n missing distinct
## 49 0 3
##
## Value fruit muesli bar yoghurt
## Frequency 33 1 15
## Proportion 0.673 0.020 0.306
## --------------------------------------------------------------------------------
## Adult.Tray.2
## n missing distinct
## 49 0 5
##
## lowest : cookies fruit muesli bar yoghurt
## highest: cookies fruit muesli bar yoghurt
##
## Value cookies fruit muesli bar yoghurt
## Frequency 17 3 9 2 18
## Proportion 0.347 0.061 0.184 0.041 0.367
## --------------------------------------------------------------------------------
## Adult.Tray.3
## n missing distinct
## 49 0 5
##
## lowest : cookies fruit muesli bar yoghurt
## highest: cookies fruit muesli bar yoghurt
##
## Value cookies fruit muesli bar yoghurt
## Frequency 17 13 4 2 13
## Proportion 0.347 0.265 0.082 0.041 0.265
## --------------------------------------------------------------------------------d <- subset(dta, select = c(id,Group,age,Q2,Q1,Q5,Q16,
Child.Tray.1:Child.Tray.3))
names(d) <- c('id','Group','age','sex','relation',
'suburb','whoBuys.food',
'child.1','child.2','child.3')
# Labels
describe(d$relation)## d$relation
## n missing distinct
## 49 0 3
##
## Value Aunty Dad Mum
## Frequency 2 22 25
## Proportion 0.041 0.449 0.510table(d$relation)##
## Aunty
## 2
## Dad
## 22
## Mum
## 25d$relation <- factor(substr(d$relation,1,1), level = c("A","D","M"),
labels = c('Aunty','Dad','Mum'))
d$Group <- factor(as.numeric(d$Group), levels = 1:3,
labels = c('Control','healthy','unhealthy'))
table(d$Group)##
## Control healthy unhealthy
## 17 16 16# classify healthy choices
d$c1 <- ifelse(d$child.1 %in% c('fruit','yogurt'), 1,0)
d$c2 <- ifelse(d$child.2 %in% c('fruit','yogurt'), 1,0)
d$c3 <- ifelse(d$child.3 %in% c('fruit','yogurt'), 1,0)
table(d$c1, d$child.1)##
## cookies fruit muesli bar yoghurt
## 0 21 0 9 9
## 1 0 10 0 0table(d$c2, d$child.2)##
##
## 0 17
## 1 0
##
## cookies fruit muesli bar yoghurt
## 0 8 0 7 1
## 1 0 16 0 0table(d$c3, d$child.3)##
##
## 0 17
## 1 0
##
## cookies fruit muesli bar
## 0 22 0 7
## 1 0 3 0head(d)## id Group age sex relation
## 1 1 unhealthy 6.362765 Female Mum
## 2 2 unhealthy 10.261465 Male Dad
## 3 3 unhealthy 10.042437 Female Mum
## 4 4 unhealthy 7.863107 Male Dad
## 5 5 unhealthy 9.738535 Female Mum
## 6 6 healthy 10.951403 Female Mum
## suburb
## 1 Narellan
## 2 Mount Annan
## 3 Spring Farm
## 4 Harrington park
## 5 Narellan
## 6 Rosemeadow
## whoBuys.food
## 1 Mum
## 2 Mum and dad
## 3 Myself
## 4 Mum and dad
## 5 Mum
## 6 Mum
## child.1 child.2 child.3 c1 c2 c3
## 1 cookies muesli bar cookies 0 0 0
## 2 cookies fruit cookies 0 1 0
## 3 cookies cookies fruit 0 0 1
## 4 cookies cookies cookies 0 0 0
## 5 muesli bar cookies cookies 0 0 0
## 6 fruit fruit fruit 1 1 1kableone(CreateTableOne(data = d,
vars = c("age","sex","relation",
"whoBuys.food"),
factorVars = c("sex","relation",
"whoBuys.food"),
strata = "Group",
addOverall = T))| Overall | Control | healthy | unhealthy | p | test | |
|---|---|---|---|---|---|---|
| n | 49 | 17 | 16 | 16 | ||
| age (mean (SD)) | 8.42 (1.97) | 8.44 (1.88) | 8.38 (1.99) | 8.44 (2.16) | 0.995 | |
| sex = Male (%) | 23 (46.9) | 9 (52.9) | 8 (50.0) | 6 (37.5) | 0.645 | |
| relation (%) | 0.349 | |||||
| Aunty | 2 ( 4.1) | 0 ( 0.0) | 0 ( 0.0) | 2 (12.5) | ||
| Dad | 22 (44.9) | 8 (47.1) | 8 (50.0) | 6 (37.5) | ||
| Mum | 25 (51.0) | 9 (52.9) | 8 (50.0) | 8 (50.0) | ||
| whoBuys.food (%) | 0.150 | |||||
| Aunty | 1 ( 2.0) | 0 ( 0.0) | 0 ( 0.0) | 1 ( 6.2) | ||
| Dad | 4 ( 8.2) | 1 ( 5.9) | 2 (12.5) | 1 ( 6.2) | ||
| Mum | 21 (42.9) | 6 (35.3) | 7 (43.8) | 8 (50.0) | ||
| Mum and dad | 15 (30.6) | 9 (52.9) | 2 (12.5) | 4 (25.0) | ||
| Myself | 1 ( 2.0) | 0 ( 0.0) | 0 ( 0.0) | 1 ( 6.2) | ||
| Parents | 7 (14.3) | 1 ( 5.9) | 5 (31.2) | 1 ( 6.2) |
# long file for repeat measure of group 2 and group 3
names(d)## [1] "id" "Group" "age" "sex" "relation"
## [6] "suburb" "whoBuys.food" "child.1" "child.2" "child.3"
## [11] "c1" "c2" "c3"d.long <- reshape(d, idvar = 'id',
varying = list(11:13),
direction = 'long',
v.names = 'yC')
d.long <- d.long[order(d.long$id, d.long$time),]
d.long$n <- 1
head(d.long)## id Group age sex relation
## 1.1 1 unhealthy 6.362765 Female Mum
## 1.2 1 unhealthy 6.362765 Female Mum
## 1.3 1 unhealthy 6.362765 Female Mum
## 2.1 2 unhealthy 10.261465 Male Dad
## 2.2 2 unhealthy 10.261465 Male Dad
## 2.3 2 unhealthy 10.261465 Male Dad
## suburb
## 1.1 Narellan
## 1.2 Narellan
## 1.3 Narellan
## 2.1 Mount Annan
## 2.2 Mount Annan
## 2.3 Mount Annan
## whoBuys.food
## 1.1 Mum
## 1.2 Mum
## 1.3 Mum
## 2.1 Mum and dad
## 2.2 Mum and dad
## 2.3 Mum and dad
## child.1 child.2 child.3 time yC n
## 1.1 cookies muesli bar cookies 1 0 1
## 1.2 cookies muesli bar cookies 2 0 1
## 1.3 cookies muesli bar cookies 3 0 1
## 2.1 cookies fruit cookies 1 0 1
## 2.2 cookies fruit cookies 2 1 1
## 2.3 cookies fruit cookies 3 0 1# Tables of outcome
# Children
stat.table(Group, contents = list(N = sum(n),
'Healthy (n)' = sum(yC),
'Healthy (%)' = ratio(yC,n, 100)),
data = d.long,
margin = T)## ------------------------------------
## Group N Healthy Healthy
## (n) (%)
## ------------------------------------
## Control 51.00 2.00 3.92
## healthy 48.00 21.00 43.75
## unhealthy 48.00 6.00 12.50
##
## Total 147.00 29.00 19.73
## ------------------------------------stat.table(sex, contents = list(N = sum(n),
'Healthy (n)' = sum(yC),
'Healthy (%)' = ratio(yC,n, 100)),
data = d.long,
margin = T)## ---------------------------------
## sex N Healthy Healthy
## (n) (%)
## ---------------------------------
## Female 78.00 16.00 20.51
## Male 69.00 13.00 18.84
##
## Total 147.00 29.00 19.73
## ---------------------------------stat.table(relation, contents = list(N = sum(n),
'Healthy (n)' = sum(yC),
'Healthy (%)' = ratio(yC,n, 100)),
data = d.long,
margin = T)## -----------------------------------
## relation N Healthy Healthy
## (n) (%)
## -----------------------------------
## Aunty 6.00 1.00 16.67
## Dad 66.00 13.00 19.70
## Mum 75.00 15.00 20.00
##
## Total 147.00 29.00 19.73
## -----------------------------------kableone(CreateTableOne(data = d.long,
vars = "yC",
factorVars = "yC",
strata = "Group",
addOverall = T))| Overall | Control | healthy | unhealthy | p | test | |
|---|---|---|---|---|---|---|
| n | 147 | 51 | 48 | 48 | ||
| yC = 1 (%) | 29 (19.7) | 2 (3.9) | 21 (43.8) | 6 (12.5) | <0.001 |
# Just Group
binGLME.adj <- glmer(cbind(yC,n) ~ Group + (1|id), family = binomial, data = d.long)
tab <- summary(binGLME.adj)
tab## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(yC, n) ~ Group + (1 | id)
## Data: d.long
##
## AIC BIC logLik deviance df.resid
## 107.3 119.3 -49.7 99.3 143
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6614 -0.3536 -0.1980 -0.1980 3.4307
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0 0
## Number of obs: 147, groups: id, 49
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.2387 0.7208 -4.493 7.03e-06 ***
## Grouphealthy 2.4120 0.7669 3.145 0.00166 **
## Groupunhealthy 1.1592 0.8409 1.379 0.16803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Grphlt
## Grouphelthy -0.940
## Groupnhlthy -0.857 0.806
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularcbind('Rate Ratio =' = round(exp(tab$coef[2:3,1]),2),
'(95% CI)' = paste(round(exp(tab$coef[2:3,1] - (1.96 * tab$coef[2:3,2])),2),',',
round(exp(tab$coef[2:3,1] + (1.96 * tab$coef[2:3,2])),2), sep = ''))## Rate Ratio = (95% CI)
## Grouphealthy "11.16" "2.48,50.15"
## Groupunhealthy "3.19" "0.61,16.57"# Group adjusted for age and sex
binGLME.adj <- glmer(cbind(yC,n) ~ Group +
age +
sex +
(1|id), family = binomial, data = d.long)
tab <- summary(binGLME.adj)
tab## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(yC, n) ~ Group + age + sex + (1 | id)
## Data: d.long
##
## AIC BIC logLik deviance df.resid
## 111.1 129.0 -49.6 99.1 141
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.7210 -0.3767 -0.2162 -0.1850 3.1330
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 3.781e-16 1.944e-08
## Number of obs: 147, groups: id, 49
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.86870 1.17215 -2.447 0.01439 *
## Grouphealthy 2.41800 0.76756 3.150 0.00163 **
## Groupunhealthy 1.14800 0.84267 1.362 0.17309
## age -0.03645 0.10849 -0.336 0.73688
## sexMale -0.14321 0.43456 -0.330 0.74174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Grphlt Grpnhl age
## Grouphelthy -0.550
## Groupnhlthy -0.525 0.804
## age -0.770 -0.033 -0.015
## sexMale -0.238 -0.014 0.055 0.088
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularcbind('Rate Ratio =' = round(exp(tab$coef[2:5,1]),2),
'(95% CI)' = paste(round(exp(tab$coef[2:5,1] - (1.96 * tab$coef[2:5,2])),2),',',
round(exp(tab$coef[2:5,1] + (1.96 * tab$coef[2:5,2])),2), sep = ''))## Rate Ratio = (95% CI)
## Grouphealthy "11.22" "2.49,50.52"
## Groupunhealthy "3.15" "0.6,16.44"
## age "0.96" "0.78,1.19"
## sexMale "0.87" "0.37,2.03"# Group, age, sex and parent
binGLME.adj <- glmer(cbind(yC,n) ~ Group +
age +
sex +
relation +
(1|id), family = binomial, data = d.long)
tab <- summary(binGLME.adj)
tab## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(yC, n) ~ Group + age + sex + relation + (1 | id)
## Data: d.long
##
## AIC BIC logLik deviance df.resid
## 114.8 138.8 -49.4 98.8 139
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.7131 -0.3744 -0.2202 -0.1915 3.0697
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.375e-16 1.172e-08
## Number of obs: 147, groups: id, 49
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6360 1.5407 -1.711 0.08709 .
## Grouphealthy 2.3590 0.7684 3.070 0.00214 **
## Groupunhealthy 1.0598 0.8644 1.226 0.22021
## age -0.0328 0.1106 -0.297 0.76677
## sexMale -15.4563 512.0000 -0.030 0.97592
## relationDad 15.1213 512.0012 0.030 0.97644
## relationMum -0.2159 1.2191 -0.177 0.85941
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Grphlt Grpnhl age sexMal rltnDd
## Grouphelthy -0.427
## Groupnhlthy -0.536 0.786
## age -0.441 -0.030 -0.052
## sexMale 0.000 0.000 0.000 0.000
## relationDad -0.001 0.000 0.000 -0.001 -1.000
## relationMum -0.649 0.010 0.221 -0.204 0.000 0.002
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularcbind('Rate Ratio =' = round(exp(tab$coef[2:7,1]),2),
'(95% CI)' = paste(round(exp(tab$coef[2:7,1] - (1.96 * tab$coef[2:7,2])),2),',',
round(exp(tab$coef[2:7,1] + (1.96 * tab$coef[2:7,2])),2), sep = ''))## Rate Ratio = (95% CI)
## Grouphealthy "10.58" "2.35,47.71"
## Groupunhealthy "2.89" "0.53,15.71"
## age "0.97" "0.78,1.2"
## sexMale "0" "0,Inf"
## relationDad "3690499.01" "0,Inf"
## relationMum "0.81" "0.07,8.79"d <- subset(dta, select = c(id,Group,age,Q2,Q1,Q5,Q16,
Adult.Tray.1:Adult.Tray.3))
names(d) <- c('id','Group','age','sex','relation',
'suburb','whoBuys.food',
'Adult.1','Adult.2','Adult.3')
# Labels
d$relation <- factor(substr(d$relation,1,1), level = c("A","D","M"),
labels = c('Aunty','Dad','Mum'))
d$Group <- factor(as.numeric(d$Group), levels = 1:3,
labels = c('Control','healthy','unhealthy'))
table(d$Group)##
## Control healthy unhealthy
## 17 16 16# classify healthy choices
d$a1 <- ifelse(d$Adult.1 %in% c('fruit','yogurt'), 1,0)
d$a2 <- ifelse(d$Adult.2 %in% c('fruit','yogurt'), 1,0)
d$a3 <- ifelse(d$Adult.3 %in% c('fruit','yogurt'), 1,0)
table(d$a1, d$Adult.1)##
## fruit muesli bar yoghurt
## 0 0 1 15
## 1 33 0 0table(d$a2, d$Adult.2)##
##
## 0 17
## 1 0
##
## cookies fruit muesli bar yoghurt
## 0 3 0 2 18
## 1 0 9 0 0table(d$a3, d$Adult.3)##
##
## 0 17
## 1 0
##
## cookies fruit muesli bar yoghurt
## 0 13 0 2 13
## 1 0 4 0 0head(d)## id Group age sex relation
## 1 1 unhealthy 6.362765 Female Mum
## 2 2 unhealthy 10.261465 Male Dad
## 3 3 unhealthy 10.042437 Female Mum
## 4 4 unhealthy 7.863107 Male Dad
## 5 5 unhealthy 9.738535 Female Mum
## 6 6 healthy 10.951403 Female Mum
## suburb
## 1 Narellan
## 2 Mount Annan
## 3 Spring Farm
## 4 Harrington park
## 5 Narellan
## 6 Rosemeadow
## whoBuys.food
## 1 Mum
## 2 Mum and dad
## 3 Myself
## 4 Mum and dad
## 5 Mum
## 6 Mum
## Adult.1 Adult.2 Adult.3 a1 a2 a3
## 1 yoghurt yoghurt cookies 0 0 0
## 2 fruit yoghurt cookies 1 0 0
## 3 fruit muesli bar cookies 1 0 0
## 4 fruit fruit muesli bar 1 1 0
## 5 yoghurt fruit cookies 0 1 0
## 6 yoghurt yoghurt yoghurt 0 0 0# long file for repeat measure of group 2 and group 3
names(d)## [1] "id" "Group" "age" "sex" "relation"
## [6] "suburb" "whoBuys.food" "Adult.1" "Adult.2" "Adult.3"
## [11] "a1" "a2" "a3"d.long <- reshape(d, idvar = 'id',
varying = list(11:13),
direction = 'long',
v.names = 'yA')
d.long <- d.long[order(d.long$id, d.long$time),]
d.long$n <- 1
describe(d$id)## d$id
## n missing distinct Info Mean Gmd .05 .10
## 49 0 49 1 25 16.67 3.4 5.8
## .25 .50 .75 .90 .95
## 13.0 25.0 37.0 44.2 46.6
##
## lowest : 1 2 3 4 5, highest: 45 46 47 48 49describe(d.long$id)## d.long$id
## n missing distinct Info Mean Gmd .05 .10
## 147 0 49 1 25 16.44 3.0 5.6
## .25 .50 .75 .90 .95
## 13.0 25.0 37.0 44.4 47.0
##
## lowest : 1 2 3 4 5, highest: 45 46 47 48 49head(d.long)## id Group age sex relation
## 1.1 1 unhealthy 6.362765 Female Mum
## 1.2 1 unhealthy 6.362765 Female Mum
## 1.3 1 unhealthy 6.362765 Female Mum
## 2.1 2 unhealthy 10.261465 Male Dad
## 2.2 2 unhealthy 10.261465 Male Dad
## 2.3 2 unhealthy 10.261465 Male Dad
## suburb
## 1.1 Narellan
## 1.2 Narellan
## 1.3 Narellan
## 2.1 Mount Annan
## 2.2 Mount Annan
## 2.3 Mount Annan
## whoBuys.food
## 1.1 Mum
## 1.2 Mum
## 1.3 Mum
## 2.1 Mum and dad
## 2.2 Mum and dad
## 2.3 Mum and dad
## Adult.1 Adult.2 Adult.3 time yA n
## 1.1 yoghurt yoghurt cookies 1 0 1
## 1.2 yoghurt yoghurt cookies 2 0 1
## 1.3 yoghurt yoghurt cookies 3 0 1
## 2.1 fruit yoghurt cookies 1 1 1
## 2.2 fruit yoghurt cookies 2 0 1
## 2.3 fruit yoghurt cookies 3 0 1# Tables of outcome
# Adults
stat.table(Group, contents = list(N = sum(n),
'Healthy (n)' = sum(yA),
'Healthy (%)' = ratio(yA,n, 100)),
data = d.long,
margin = T)## ------------------------------------
## Group N Healthy Healthy
## (n) (%)
## ------------------------------------
## Control 51.00 8.00 15.69
## healthy 48.00 18.00 37.50
## unhealthy 48.00 20.00 41.67
##
## Total 147.00 46.00 31.29
## ------------------------------------stat.table(sex, contents = list(N = sum(n),
'Healthy (n)' = sum(yA),
'Healthy (%)' = ratio(yA,n, 100)),
data = d.long,
margin = T)## ---------------------------------
## sex N Healthy Healthy
## (n) (%)
## ---------------------------------
## Female 78.00 24.00 30.77
## Male 69.00 22.00 31.88
##
## Total 147.00 46.00 31.29
## ---------------------------------stat.table(relation, contents = list(N = sum(n),
'Healthy (n)' = sum(yA),
'Healthy (%)' = ratio(yA,n, 100)),
data = d.long,
margin = T)## -----------------------------------
## relation N Healthy Healthy
## (n) (%)
## -----------------------------------
## Aunty 6.00 2.00 33.33
## Dad 66.00 22.00 33.33
## Mum 75.00 22.00 29.33
##
## Total 147.00 46.00 31.29
## -----------------------------------kableone(CreateTableOne(data = d.long,
vars = "yA",
factorVars = "yA",
strata = "Group",
addOverall = T))| Overall | Control | healthy | unhealthy | p | test | |
|---|---|---|---|---|---|---|
| n | 147 | 51 | 48 | 48 | ||
| yA = 1 (%) | 46 (31.3) | 8 (15.7) | 18 (37.5) | 20 (41.7) | 0.011 |
# Just Group
binGLME.adj <- glmer(cbind(yA,n) ~ Group + (1|id), family = binomial, data = d.long)
tab <- summary(binGLME.adj)
tab## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(yA, n) ~ Group + (1 | id)
## Data: d.long
##
## AIC BIC logLik deviance df.resid
## 150.8 162.8 -71.4 142.8 143
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6455 -0.6124 -0.3961 0.6390 1.5053
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0 0
## Number of obs: 147, groups: id, 49
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.8524 0.3803 -4.871 1.11e-06 ***
## Grouphealthy 0.8716 0.4701 1.854 0.0637 .
## Groupunhealthy 0.9769 0.4642 2.105 0.0353 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Grphlt
## Grouphelthy -0.809
## Groupnhlthy -0.819 0.663
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularcbind('Rate Ratio =' = round(exp(tab$coef[2:3,1]),2),
'(95% CI)' = paste(round(exp(tab$coef[2:3,1] - (1.96 * tab$coef[2:3,2])),2),',',
round(exp(tab$coef[2:3,1] + (1.96 * tab$coef[2:3,2])),2), sep = ''))## Rate Ratio = (95% CI)
## Grouphealthy "2.39" "0.95,6.01"
## Groupunhealthy "2.66" "1.07,6.6"# Group adjusted for age and sex
binGLME.adj <- glmer(cbind(yA,n) ~ Group +
age +
sex +
(1|id), family = binomial, data = d.long)
tab <- summary(binGLME.adj)
tab## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(yA, n) ~ Group + age + sex + (1 | id)
## Data: d.long
##
## AIC BIC logLik deviance df.resid
## 154.6 172.6 -71.3 142.6 141
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6899 -0.6117 -0.4048 0.6243 1.6470
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0 0
## Number of obs: 147, groups: id, 49
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.13456 0.85252 -2.504 0.0123 *
## Grouphealthy 0.87903 0.47081 1.867 0.0619 .
## Groupunhealthy 0.98122 0.46689 2.102 0.0356 *
## age 0.02783 0.08931 0.312 0.7554
## sexMale 0.08802 0.34623 0.254 0.7993
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Grphlt Grpnhl age
## Grouphelthy -0.390
## Groupnhlthy -0.342 0.661
## age -0.867 0.023 -0.049
## sexMale -0.135 0.040 0.096 -0.099
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularcbind('Rate Ratio =' = round(exp(tab$coef[2:5,1]),2),
'(95% CI)' = paste(round(exp(tab$coef[2:5,1] - (1.96 * tab$coef[2:5,2])),2),',',
round(exp(tab$coef[2:5,1] + (1.96 * tab$coef[2:5,2])),2), sep = ''))## Rate Ratio = (95% CI)
## Grouphealthy "2.41" "0.96,6.06"
## Groupunhealthy "2.67" "1.07,6.66"
## age "1.03" "0.86,1.22"
## sexMale "1.09" "0.55,2.15"# Group, age, sex and parent
binGLME.adj <- glmer(cbind(yA,n) ~ Group +
age +
sex +
relation +
(1|id), family = binomial, data = d.long)
tab <- summary(binGLME.adj)
tab## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(yA, n) ~ Group + age + sex + relation + (1 | id)
## Data: d.long
##
## AIC BIC logLik deviance df.resid
## 157.7 181.6 -70.8 141.7 139
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6924 -0.6142 -0.4189 0.6034 1.5673
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0 0
## Number of obs: 147, groups: id, 49
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.18454 1.10783 -1.972 0.0486 *
## Grouphealthy 0.81639 0.47193 1.730 0.0837 .
## Groupunhealthy 0.93896 0.47664 1.970 0.0488 *
## age 0.02167 0.09098 0.238 0.8118
## sexMale -15.70910 1448.15502 -0.011 0.9913
## relationDad 15.96804 1448.15556 0.011 0.9912
## relationMum 0.15725 0.89893 0.175 0.8611
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Grphlt Grpnhl age sexMal rltnDd
## Grouphelthy -0.296
## Groupnhlthy -0.384 0.647
## age -0.522 0.032 -0.082
## sexMale 0.000 0.000 0.000 0.000
## relationDad 0.000 0.000 0.000 0.000 -1.000
## relationMum -0.640 -0.011 0.196 -0.204 0.000 0.001
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularcbind('Rate Ratio =' = round(exp(tab$coef[2:7,1]),2),
'(95% CI)' = paste(round(exp(tab$coef[2:7,1] - (1.96 * tab$coef[2:7,2])),2),',',
round(exp(tab$coef[2:7,1] + (1.96 * tab$coef[2:7,2])),2), sep = ''))## Rate Ratio = (95% CI)
## Grouphealthy "2.26" "0.9,5.71"
## Groupunhealthy "2.56" "1,6.51"
## age "1.02" "0.85,1.22"
## sexMale "0" "0,Inf"
## relationDad "8606616.74" "0,Inf"
## relationMum "1.17" "0.2,6.82"#-------------------------------- The End -------------------------------------------------------