First we need to set our working directory and load our data into R

setwd("C:/Users/jedia/Desktop/RWD")
USstates<-read.csv("C:/Users/jedia/Desktop/RWD/USstates.csv")

Next we’re going to create objects for our key variables, for ease of access

urban<-USstates$urban
secular<-USstates$secularism
relig<-USstates$relig_high

Now we’re going to make simple scatterplots to be sure there’s some kind of correlation between our variables

plot(urban~secular)

plot(urban~relig)

It looks like there’s some kind of relationship, so we’ll use a function to actually compare these and show us if there’s any significance

urbsecular<-lm(urban~secular, USstates)
urbrelig<-lm(urban~relig, USstates)
summary(urbsecular)
## 
## Call:
## lm(formula = urban ~ secular, data = USstates)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.131  -9.330   2.077  11.172  21.629 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 63.83297    4.32024  14.775   <2e-16 ***
## secular      0.08054    0.03902   2.064   0.0444 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.43 on 48 degrees of freedom
## Multiple R-squared:  0.08155,    Adjusted R-squared:  0.06241 
## F-statistic: 4.262 on 1 and 48 DF,  p-value: 0.04441
summary(urbrelig)
## 
## Call:
## lm(formula = urban ~ relig, data = USstates)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -41.127  -9.123   1.884  12.793  22.750 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  86.5093     9.4221   9.182 3.84e-12 ***
## relig        -0.3761     0.2333  -1.612    0.114    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.67 on 48 degrees of freedom
## Multiple R-squared:  0.05135,    Adjusted R-squared:  0.03159 
## F-statistic: 2.598 on 1 and 48 DF,  p-value: 0.1135

It looks as though there is a slight significant relationship present between urban populations and secularism, but there is no significant relationship between urban populations and high religiosity. Nevertheless, we’re going to make figures for both.

The first figure we need to make is a table showing our information. For this, we’re going to use an r package called gt summary. install.packages(“gtsummary”)

library(gtsummary)
bold_p(x=modify_header(x=tbl_regression(urbsecular, intercept=TRUE, label=list(secular ~ "Secularism")), update = list(estimate ~ "**Coefficient**")))
Characteristic Coefficient 95% CI1 p-value
(Intercept) 64 55, 73 <0.001
Secularism 0.08 0.00, 0.16 0.044

1 CI = Confidence Interval 

bold_p(x=modify_header(x=tbl_regression(urbrelig, intercept=TRUE, label=list(relig ~ "Religiosity")), update = list(estimate ~ "**Coefficent**")))
Characteristic Coefficent 95% CI1 p-value
(Intercept) 87 68, 105 <0.001
Religiosity -0.38 -0.85, 0.09 0.11

1 CI = Confidence Interval

Now that we’ve made those, its time to make figures for our graphs.

We’re first going to use the predict(), so that we can apply our confidence intervals to our actual graphs.

secupredict<-predict(urbsecular, interval="confidence")
religpredict<-predict(urbrelig, interval="confidence")

These will be fundamental to creating our confidence interval representation, as they provide y values for every potential x value, and will be the borders of our polygon function.

As the last step before we make our figures, we need to create sorted versions of our variables to use for the lines and polygon

secu<-sort(secular,index.return=T)$ix
reli<-sort(relig,index.return=T)$ix

Now we just need to create our scatterplots, add a line of best fit using our regression model, and add the confidence intervals using our prediction model

plot(urban~secular,col="dark Red", main="Secularism by Urban population", xlab="Secularism Scale", ylab="% Urban Population", pch=16)
mtext("p=.0440", side=3)
lines(secular[secu], secupredict[secu, 1], col = "Red", lwd=2)
polygon(c(rev(secular[secu]), secular[secu]), c(rev(secupredict[secu,3]), secupredict[secu, 2]), col = rgb(.7,.7,.7,.4), border=NA)

plot(urban~relig,col="Dark Blue", main="Religiosity by Urban Population", xlab= "% High Religiosity", ylab="% Urban Population", pch=16)
mtext("p=.114", side=3)
lines(relig[reli], religpredict[reli, 1], col= "Light Blue", lwd=2)
polygon(c(rev(relig[reli]), relig[reli]), c(rev(religpredict[reli,3]), religpredict[reli, 2]), col = rgb(.7,.7,.7,.4), border= NA)

Now we just need to export our figures and upload them to our file