enquirysugar911.Rmd

---
title: " Blood Pressure in NDNS"
output:
  html_notebook: default
  word_document:
    toc: yes
    fig_caption: yes
  html_document:
    toc: yes
    df_print: paged
  pdf_document: default
---


```{r library setup, message=FALSE, warning=FALSE, include=FALSE}
library(haven)
library(data.table)
library(ggplot2)

sav1 <- "UKDA-6533-spss/spss/spss25/ndns_rp_yr9-11a_personleveldietarydata_uk_20210831.sav"
sav2 <- "UKDA-6533-spss/spss/spss25/ndns_rp_yr9-11a_indiv_20211020.sav"
#sav3 <- "UKDA-6533-spss/spss/spss25/ndns_yr1-3indiva_weights.sav"
```

####   What is the best model of Systolic and Diastolic Blood pressure in the NDNs dataset? And How should this influence policy?


In a representative UK population (NDNS) what is the best model for understanding the relationship between diet and health, between intake of fundamental nutrients and blood pressure?

## Introduction 

Sodium is a key electrolyte in cellular physiology. One core function of the cell is to maintain a sodium concentration gradient across the cell wall.  Maintaining sodium levels is therefore an essential part of all higher organisms. This role is largely taken by the kidney and modified by receptors and hormones from across the organism. It also relies on ingestion and taste, or diet. 

I would like to understand what to tell my patients. Current advice is to reduce sodium intake but it can be difficult to identify the best way to do this. Recent critiques of the role of ingested sodium in blood pressure have looked again at the best form for this advice.  The National dietary nutrition dataset is a rolling cross sectional study with linked data on ingestion of sodium and of blood pressure. It also contains data on potential confounding factors such as age sex weight race and income.

Machine learning is a way of interrogating data sets to identify potential models to explain and predict an outcome. In this case the outcome is systolic or diastolic blood pressure.

Causation analysis working with regression models can help to better identify the directional role of particular variables in a model. This helps to identify which variables to include in models and which combination is most likely to give a clinically significant answer.



```{r read data, message=FALSE, warning=FALSE, include=FALSE}
#savx <- "UKDA-6533-spss/spss/spss25/"

sav1d <- read_sav(sav1)
#View(sav1d)

sav2d <- read_sav(sav2)
View(sav2d)

#savxd <- read_sav(savx)
 
#sav3d <- read_sav(sav3)



#yr1-4a person level dietary data_uk_v2 data
persdat <- subset(sav1d[c("seriali","AgeR", "Sex","Country","Sodiummg","Calciummg","TotalEMJ")])
sugarset <- subset(sav1d[c("seriali", "Totalsugarsg","Glucoseg","Sucroseg","Fructoseg","Lactoseg","SOFTDRINKSLOWCALORIE","SOFTDRINKSNOTLOWCALORIE", "TEACOFFEEANDWATER")])

persdat <- as.data.table(persdat)
sugarset <- as.data.table(sugarset)

persdat[, Sex := factor(Sex, levels = 1:2, labels = c("Male", "Female"))]
persdat[, Country := factor(Country)]

#summary(persdat)
#summary(sugarset)
```




```{r setup variables, message=FALSE, warning=FALSE, include=FALSE}
 #measured and recorded data yr1-4a_indiv_uk 
#subsets of the table to identify grouped information
bpset <- subset(sav2d[c("seriali","omsysval","omdiaval")])
ethnset <- subset(sav2d[c("seriali" ,"ethgrp5","ethgrp2")])
saltset <- subset(sav2d[c("seriali","SaltChk","SalHowC","SltSHow" )])
#medsset <- subset(sav2d[c("seriali", "bpmedc","bpmedd")])
hypset <- subset(sav2d[c("seriali", "hyper140_2", "hibp140_2", "hyper1_2", "highbp1_2")])
incset <- subset(sav2d[c("seriali","nssec8")])
measset <- subset(sav2d[c("seriali","htval2","wtval2","bmival2","bmivg5")])
ageset <- subset(sav2d[c("seriali","agegad1", "agegad2", "agegch1","agegr1")])
```



```{r rearrange, message=FALSE, warning=FALSE, include=FALSE}
#change subsets to datatables  
 bpset <- as.data.table(bpset)
 saltset <- as.data.table(saltset)
 #medsset <- as.data.table(medsset)
 hypset <- as.data.table(hypset)
 ethnset <- as.data.table(ethnset)
 incset <- as.data.table(incset)
 measset <- as.data.table(measset)
 ageset <- as.data.table(ageset)
 
 #define factors 
 #saltset
 saltset[, SaltChk := factor(SaltChk, levels = 1:8, labels = c("Salt",
	"Salt substitute",
	"Neither",
	"Item not applicable",
	"No answer/refused",
	"Don't know",
	"Qn not applicable to survey year",
	"Schedule not applicable"))]
 saltset[, SalHowC := factor(SalHowC, levels = 1:8,labels = c("Always",
"Usually",
"Sometimes",
	",Item not applicable",
	"No answer/refused",
	" Don't know",
"Qn not applicable to survey year",
"Schedule not applicable"))]
 saltset[, SltSHow := factor(SltSHow,levels = 1:8,labels = c("Always",
                                                            "Usually",
                                                            "Sometimes",
                                                            ",Item not applicable",
                                                            "No answer/refused",
                                                            " Don't know",
                                                            "Qn not applicable to survey year",
                                                            "Schedule not applicable"))]
 
#ethnset
ethnset[ , ethgrp5 := factor(ethgrp5, levels = 1:5, labels = c( 'White'
                                                             , 'Mixed ethnic group'
                                                             , 'Black or Black British'
                                                             , 'Asian or asian British'
                                                             , 'Any other group'))]
ethnset[ , ethgrp2 := factor(ethgrp2, levels = 1:2, labels = c( 'White'
                                                             , 'Non-white'))]

 
 
#bpset 
 
 bpset[, "omsysval" := as.numeric(omsysval)]
 #bpset[,Sys := as.numeric(Sys)]

 

#incset
incset[ , nssec8 := factor(nssec8, levels = 1:9, labels = c(  "Higher managerial and professional occupations", "Lower managerial and professional occupations"
               , "Intermediate occupations"
                           , "Small employers and own account workers"
                           , "Lower supervisory and technical occupations"
                           , "Semi-routine occupations"
                           , "Routine occupations"
                           , "Never worked"
                           , "Other"))]


 #hypeset
 hypset[,hyper140_2 := factor(hyper140_2,levels = 1:9, labels = c(" Normotensive untreated",
" Normotensive treated",
	"Hypertensive treated",
"Hypertensive untreated",
"No answer/refused",
	"Don't know",
"Refused, attempted but not obtained, not attempted",
	 "Qn not applicable to survey year",
"Item not applicable")) ]
 hypset[, hibp140_2 := factor(hibp140_2, levels = 1:7, labels = c("Not high BP",
"High BP",
	"No answer/refused",
	"Don't know",
	"Refused, attempted but not obtained, not attempted",
	"Qn not applicable to survey year",
	"Item not applicable")) ]
 hypset[, hyper1_2 := factor(hyper1_2,levels = 1:9, labels = c(" Normotensive untreated",
                                                          " Normotensive treated",
                                                          "Hypertensive treated",
                                                          "Hypertensive untreated",
                                                          "No answer/refused",
                                                          "Don't know",
                                                          "Refused, attempted but not obtained, not attempted",
                                                          "Qn not applicable to survey year",
                                                          "Item not applicable"))]
 hypset[, highbp1_2 := factor(highbp1_2, levels = 1:7, labels = c("Not high BP",
                                                              "High BP",
                                                              "No answer/refused",
                                                              "Don't know",
                                                              "Refused, attempted but not obtained, not attempted",
                                                              "Qn not applicable to survey year",
                                                              "Item not applicable")) ]

 
#measset
 
 measset[,bmivg5 := factor(bmivg5, levels = 1:10, labels = c("Under 18.5", "18.5 and below 25","25 and below 30","30 and below 40","Over 40","Refused","Don't know","Different variable this survey year","Not applicable to survey year","Not applicable"))]



 
```

# literature

The importance of sodium as the extracellular electrolyte was established by physiologists such as Guyton. They also demonstrated that sodium chloride infusion affected renal blood flow and blood pressure. 

Pharmacologists working on hypertension have used diuretic drugs to reduce blood pressure. These increase renal excretion. Loop diuretics work on sodium in the loop of Henle. Spironolactone works to block the effect of the natural hormone aldosterone. Aldosterone is the active member of the renin- angiotensin hormone system which plays an important role in blood pressure regulation. The most widely used treatments for hypertension work on this system blocking the angiotensin converting enzyme.

It is logical therefore that the role of sodium is important in managing blood pressure in everybody.

National governments and WHO have recommended keeping oral intake below 2.3g. However daily sodium intake seems to have been increasing. 

The NDNS time trend analysis showed changes in salt intake. Whilst BP measurements were made the results are not reported on in the paper. 

It would be useful to understand if the predicted improvements in BP were found in the same population. 

Then the meaning of these results needs to be understandable at the level of those able to implement policy. The NDNS population is structured to match the UK age and sex profile.  





```{r data summary, echo=FALSE, message=FALSE, warning=FALSE}
summary(persdat)
summary(bpset)
summary(saltset)
summary(medsset)
summary(hypset)
summary(ageset)
summary(sugarset)

#combine data from tables

persugar <- merge(persdat,sugarset, by = "seriali")
persethnsugar <- merge(ethnset , persugar, by = "seriali")
incethsugar <- merge(persethnsugar, incset , by = "seriali")


meaag <- merge(measset,ageset, by="seriali")
meaage <- merge(saltset,meaag, by ="seriali" )
meaages <- merge(meaage,hypset, by = "seriali")


nearlyalldata <- merge(incethsugar, bpset, by ="seriali")
alldata <- merge(nearlyalldata, meaages, by ="seriali")

alldatam <- alldata[Sex == "Male"]
alldataf <- alldata[Sex == "Female"]



View(alldata)

```

#method
The NDNS is a postcode randomised survey which approaches approximately 1000 people each year. 
The sample selects 500 adults and 500 children. The numbers are managed to deliver a representative sample for the UK. 

These participants are asked some basic questions and if they agree to take part are given a 4 day dietary diary. The recorded intake is then analysed for intake by food. These foods have defined contents which can then be reduced to their elemental constituents.


# the data

```{r hist1, echo=FALSE, message=FALSE, warning=FALSE}

#description of dataset 

 

hist(alldata[,omsysval])
alldata[, hist(omsysval, prob = TRUE)] # histogram
alldata[, curve(
  dnorm(x, mean(omsysval, na.rm = TRUE), sd(omsysval, na.rm = TRUE)),
  add = TRUE)] # superimpose a Normal distribution

hist(alldata[,omdiaval])



```






```{r hist2, echo=FALSE, message=FALSE, warning=FALSE}
hist(alldata[,Sodiummg])
hist(alldata[,Calciummg])
hist(alldata[,TotalEMJ])
hist(alldata[,Totalsugarsg])
hist(alldata[,wtval2])

hist(alldata[,bmival2])


dietnagraph <- ggplot(alldata, aes(Sex,Sodiummg))+ geom_boxplot()
dietnagraph

dietcagraph <- ggplot(alldata, aes(Sex,Calciummg))+ geom_boxplot()
dietcagraph

dietSugraph <- ggplot(alldata, aes(Sex,Totalsugarsg))+ geom_boxplot()
dietSugraph

dietJgraph <- ggplot(alldata, aes(Sex,TotalEMJ))+ geom_boxplot()
dietJgraph

dietJbmigraph <- ggplot(alldata, aes(bmivg5,TotalEMJ))+ geom_boxplot()
dietJbmigraph

dietJethgraph <- ggplot(alldata, aes(nssec8,TotalEMJ))+ geom_boxplot()
dietJethgraph

```








```{r graphs1, echo=FALSE, message=FALSE, warning=FALSE}


# graph section view the data
#bpdietdat[Sex == "2" ,plot("Sodiummg","Sys")]
graph1 <- ggplot(alldata, aes(Sodiummg, omsysval, colour = factor(Country)) ) + geom_smooth()
graph1
graph1a <- ggplot(alldata, aes(Sodiummg, omsysval, colour = factor(Sex)) )+ geom_smooth()
graph1a

graph1b <- ggplot(alldata, aes(Sodiummg, omsysval, colour = factor(agegad1)) ) + geom_smooth()
graph1b

graph1bm <- ggplot(alldatam, aes(Sodiummg, omsysval, colour = factor(agegad1)) ) + geom_smooth()
graph1bm

graph1bf <- ggplot(alldataf, aes(Sodiummg, omsysval, colour = factor(agegad1)) ) + geom_smooth()
graph1bf


graph1bd <- ggplot(alldata, aes(Sodiummg, omdiaval, colour = factor(agegad1)) ) + geom_smooth()
graph1bd

graph1c <- ggplot(alldata, aes(Sodiummg, omsysval, colour = factor(agegad2)) ) + geom_smooth()
graph1c

graph1bm <- ggplot(alldatam, aes(Sodiummg, omsysval, colour = factor(agegad2)) ) + geom_smooth()
graph1bm

graph1bf <- ggplot(alldataf, aes(Sodiummg, omsysval, colour = factor(agegad2)) ) + geom_smooth()
graph1bf

graph1cd <- ggplot(alldata, aes(Sodiummg, omdiaval, colour = factor(agegad2)) ) + geom_smooth()
graph1cd



```

# results

```{r graphs2, echo=FALSE, message=FALSE, warning=FALSE}
#bpdietdatq <- subset(bpdietdat[,c("seriali","Age","Sodiummg","eqvinc","omsysval", "Dias", "omsysval2", "Dias2","omsysval","omdiaval")])
#bpdietdatq <- bpdietdatq[,is.na(bpdietdat[Sys]) ]
#View(bpdietdatq)
#Summary(bpdietdatq)

alldata[, hist(wtval2, prob = TRUE)] # histogram
alldata[, curve(
  dnorm(x, mean(wtval2, na.rm = TRUE), sd(wtval2, na.rm = TRUE)),
  add = TRUE
)] # superimpose a Normal distribution

# graph section view the data
#bpdietdat[Sex == "2" ,plot("Sodiummg","Sys")]
graphw1 <- ggplot(alldata, aes(wtval2, omsysval, colour = factor(Country)) ) + geom_smooth()
graphw1
graphw1a <- ggplot(alldata, aes(wtval2, omsysval, colour = factor(Sex)) )+ geom_smooth()
graphw1a

graph1wb <- ggplot(alldata, aes(wtval2, omsysval, colour = factor(agegad1)) ) + geom_smooth()
graph1wb

graphw1bm <- ggplot(alldatam, aes(wtval2, omsysval, colour = factor(agegad1)) ) + geom_smooth()
graphw1bm

graphw1bf <- ggplot(alldataf, aes(wtval2, omsysval, colour = factor(agegad1)) ) + geom_smooth()
graphw1bf


graphw1bd <- ggplot(alldata, aes(wtval2, omdiaval, colour = factor(agegad1)) ) + geom_smooth()
graphw1bd

graphw1c <- ggplot(alldata, aes(wtval2, omsysval, colour = factor(agegad2)) ) + geom_smooth()
graphw1c

graphw1bm <- ggplot(alldatam, aes(wtval2, omsysval, colour = factor(agegad2)) ) + geom_smooth()
graphw1bm

graphw1bf <- ggplot(alldataf, aes(wtval2, omsysval, colour = factor(agegad2)) ) + geom_smooth()
graphw1bf

graphw1cd <- ggplot(alldata, aes(wtval2, omdiaval, colour = factor(agegad2)) ) + geom_smooth()
graphw1cd
```



```{r}
# graph section view the data
#bpdietdat[Sex == "2" ,plot("Sodiummg","Sys")]
graph3 <- ggplot(alldata, aes(TotalEMJ, bmival2, colour = factor(Country)) ) + geom_smooth()
graph3
graph3a <- ggplot(alldata, aes(TotalEMJ, bmival2, colour = factor(Sex)) )+ geom_smooth()
graph3a

graph3b <- ggplot(alldata, aes(TotalEMJ, bmival2, colour = factor(agegad1)) ) + geom_smooth()
graph3b

graph3bm <- ggplot(alldatam, aes(TotalEMJ, bmival2, colour = factor(agegad1)) ) + geom_smooth()
graph3bm

graph3bf <- ggplot(alldataf, aes(TotalEMJ, bmival2, colour = factor(agegad1)) ) + geom_smooth()
graph3bf


graph3bd <- ggplot(alldata, aes(TotalEMJ, bmival2, colour = factor(agegad1)) ) + geom_smooth()
graph3bd

graph3c <- ggplot(alldata, aes(TotalEMJ, bmival2, colour = factor(agegad2)) ) + geom_smooth()
graph3c

graph3bm <- ggplot(alldatam, aes(TotalEMJ, bmival2,colour = factor(agegad2)) ) + geom_smooth()
graph3bm

graph3bf <- ggplot(alldataf, aes(TotalEMJ, bmival2, colour = factor(agegad2)) ) + geom_smooth()
graph3bf

graph3cd <- ggplot(alldata, aes(TotalEMJ, bmival2, colour = factor(agegad2)) ) + geom_smooth()
graph3cd
```




```{r graphs3, echo=FALSE, message=FALSE, warning=FALSE}

graph2 <- ggplot(alldatam, aes(AgeR, Sodiummg, colour = factor(Sex)) )+ geom_smooth()
graph2

graph2a <- ggplot(alldatam, aes(AgeR, omsysval) ) + geom_smooth()
graph2a

graph2b <- ggplot(alldatam, aes(AgeR, omdiaval) ) + geom_smooth()
graph2b

#graph3 <- ggplot(alldatam, aes(Sodiummg, eqvinc) )+ geom_smooth()
#graph3

#graph4 <- ggplot(alldatam, aes(omsysval, eqvinc) ) +geom_smooth()
#graph4



#show persugstat

graph5 <- ggplot(alldatam, aes(AgeR, Fructoseg)) + geom_smooth()
graph5


#show ethnsugar

graph6 <- ggplot(alldatam, aes (AgeR , Fructoseg, colour = ethgrp2))+geom_col()
graph6

```



```{r statspart1, echo=FALSE, message=FALSE, warning=FALSE}
# statistical analysis section 

# comparison tables

incethtbl <- alldata[,table(nssec8, ethgrp2)]
incethtbl

incethtbl2 <- alldata[,table(nssec8, ethgrp5)]
incethtbl2

agetbl1 <-alldata[,table(agegr1,Sex)]
agetbl1




#how many are on bp medication ? 
#medstbl1 <- alldata[ , table(bpmedc,bpmedd)]
#medstbl1

#hyptabl1 <- alldata[,table(hibp140_2, bpmedc)]
#hyptabl1

#hyptabl2 <- alldata[,table(hyper1_2, bpmedc)]
#hyptabl2

#hyptabl3 <- alldata[,table(hyper140_2, bpmedc)]
#hyptabl3

#hyptabl4 <- alldata[,table(highbp1_2, bpmedc)]
#hyptabl4
```

# analysis


```{r analysis1, echo=FALSE}
#random correlations

meds <- t.test(medstbl1 ,  na.rm = TRUE)

incSys <-  t.test( hyptabl1,  na.rm = TRUE)
incSys
meds




#linear regression models
lm1 <- lm(omsysval ~ agegad2 + Sex +Sodiummg + Calciummg + Fructoseg  + ethgrp2 +TotalEMJ +wtval2 , alldata)
summary(lm1)
plot(lm1)

lm1m <- lm(omsysval ~ agegad2  +Sodiummg + Calciummg + Fructoseg  + ethgrp2 +TotalEMJ +wtval2 , alldatam)
summary(lm1m)
plot(lm1m)

lm1f <- lm(omsysval ~ agegad2  +Sodiummg + Calciummg + Fructoseg  + ethgrp2 +TotalEMJ +wtval2, alldataf)
summary(lm1f)
plot(lm1f)


lmD1 <- lm(omdiaval ~ agegad2 + Sex +Sodiummg + Calciummg + Fructoseg  + ethgrp2 +TotalEMJ +wtval2,  alldata)
summary(lmD1)

lm2 <- lm(omsysval ~  Sex + AgeR + TEACOFFEEANDWATER, alldata)
summary(lm2)
plot(lm2)

lmD2 <- lm(omdiaval ~  Sex + AgeR + TEACOFFEEANDWATER, alldata)
summary(lmD2)
```




```{r mach1, echo=TRUE}
library(caret)
set.seed(42)

alldataSys <- alldata[!is.na(alldata$omsysval),]
alldataSys <- as.data.frame(alldataSys)
simplealldata <- subset(alldataSys[c("omsysval", "AgeR","wtval2","Sex")])

                      
modelSys <- train( omsysval~ ., simplealldata,
               method = "lm",
               trControl = trainControl(method = "cv",
                                        number = 10,
                                        verboseIter = TRUE), na.action = na.pass)
modelSys

```

```{r mach2, echo=TRUE}
p <- predict(modelSys, simplealldata)
error <- p - simplealldata$omsysval
rmse_xval <- sqrt(mean(error^2)) ## xval RMSE
rmse_xval
```

#discussion

The dataset shows that sodium intake is related to systolic blood pressure in particular age groups of men. There is much less of a relationship in women. 

These relationships are different in different age groups. 

The data support consideration of a more complex approach to preventing blood pressure. 
The simplest message for avoidance is to loose weight. Changing age and Sex are more complex. 
Salt reduction is an issue for men between 16 and 50, but even here it is more important to maintain a healthy weight than to reduce sodium intake. 
This may be due to the bodies ability to auto regulate sodium to much higher levels than oral intake allows.