change in perceived wellbeing

library(visreg)
library(broom)
library(car)

Loading required package: carData

setwd("C:/Users/Chia chen/Desktop/Brenda")

data<-read.csv("BrisSydforR.csv")

data$City<-as.factor(data$City)
summary(as.factor(data$HealthChangeYear))

  -2   -1    0    1    2 
  47  312 1313  347   65 

data$HealthChangeYearGroup<-ifelse(data$HealthChangeYear <0, 'Negative', ifelse(data$HealthChangeYear > 0, 'Positive', 'Same'))
data$HealthChangeYearGroup<-as.factor(data$HealthChangeYearGroup)
data$HealthChangeYearF<-ifelse(data$HealthChangeYear <0, '-1', ifelse(data$HealthChangeYear > 0, '1', '0'))
data$HealthChangeYearF<-as.numeric(data$HealthChangeYearF)

data$TimeWeekMoreLessF<-as.factor(data$TimeWeekMoreLess)

data$X.10minYardWeek[is.na(data$X.10minYardWeek)]<-0
data$TimeYardWeek[is.na(data$TimeYardWeek)]<-0

#hourslastweek less than 0 is NA, and more than 168 is NA 
data$HoursLastWeek<-ifelse(data$HoursLastWeek <0|data$HoursLastWeek >168, NA, data$HoursLastWeek)
#number of visit in a week is 0, then hours spent should be 0
data$HoursLastWeek<-ifelse(data$NoTimeVisitGSWeek ==0,0, data$HoursLastWeek)
data$NoTimeVisitGSWeek<-ifelse(data$HoursLastWeek ==0,0, data$NoTimeVisitGSWeek)

length(data[data$HoursLastWeek>9.99,])

[1] 50

#summary(data$HoursLastWeek)

data$HoursLastWeekGroup<-NA
#data$HoursLastWeekGroup[data$HoursLastWeek>6.01]<-7
data$HoursLastWeekGroup[data$HoursLastWeek>9.99]<-10#10 or more 
data$HoursLastWeekGroup[data$HoursLastWeek<10]<-9#9 hrs
data$HoursLastWeekGroup[data$HoursLastWeek<9]<-8#8hr
data$HoursLastWeekGroup[data$HoursLastWeek<8]<-7#7hr
data$HoursLastWeekGroup[data$HoursLastWeek<7]<-6#6hr
data$HoursLastWeekGroup[data$HoursLastWeek<6]<-5#5hr
data$HoursLastWeekGroup[data$HoursLastWeek<5]<-4#4hr
data$HoursLastWeekGroup[data$HoursLastWeek<4]<-3#3hr
data$HoursLastWeekGroup[data$HoursLastWeek<3]<-2#2hr
data$HoursLastWeekGroup[data$HoursLastWeek<2]<-1#1hr
data$HoursLastWeekGroup[data$HoursLastWeek== 0]<-0

#summary(data)

Positive change in wellbeing

positive<-data[!data$HealthChangeYearGroup=="Negative",]
#the higher number means more likely to have positive impact
m1<-glm(HealthChangeYearF~TimeYardYear+TimeGSYear+
          TimeYardWeek+X.10minYardWeek+FreqVisitGS+AvgNR+HoursLastWeekGroup+
          Age+Gender+Income+City,data=positive,family=binomial)
summary(m1)

Call:
glm(formula = HealthChangeYearF ~ TimeYardYear + TimeGSYear + 
    TimeYardWeek + X.10minYardWeek + FreqVisitGS + AvgNR + HoursLastWeekGroup + 
    Age + Gender + Income + City, family = binomial, data = positive)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8821  -0.7414  -0.5518  -0.2156   2.5907  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)        -3.47054    0.53818  -6.449 1.13e-10 ***
TimeYardYear        0.02150    0.08631   0.249   0.8033    
TimeGSYear          0.54813    0.09896   5.539 3.05e-08 ***
TimeYardWeek       -0.09831    0.04945  -1.988   0.0468 *  
X.10minYardWeek     0.13752    0.05861   2.346   0.0190 *  
FreqVisitGS         0.04380    0.03866   1.133   0.2573    
AvgNR               0.64240    0.13056   4.920 8.63e-07 ***
HoursLastWeekGroup  0.01829    0.02165   0.845   0.3984    
Age                -0.17155    0.02293  -7.482 7.34e-14 ***
Gender              0.20895    0.14320   1.459   0.1445    
Income              0.02965    0.02663   1.113   0.2655    
CitySydney         -0.11759    0.13945  -0.843   0.3991    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1462.2  on 1315  degrees of freedom
Residual deviance: 1304.4  on 1304  degrees of freedom
  (409 observations deleted due to missingness)
AIC: 1328.4

Number of Fisher Scoring iterations: 4

Negative change in wellbeing

negative<-data[!data$HealthChangeYearGroup=="Positive",]
#the higher number means more likely to have negative impact
m2<-glm((HealthChangeYearF*-1)~TimeYardYear+TimeGSYear+
          TimeYardWeek+X.10minYardWeek+FreqVisitGS+AvgNR+HoursLastWeekGroup+
          Age+Gender+Income+City,data=negative,family=binomial)
summary(m2)

Call:
glm(formula = (HealthChangeYearF * -1) ~ TimeYardYear + TimeGSYear + 
    TimeYardWeek + X.10minYardWeek + FreqVisitGS + AvgNR + HoursLastWeekGroup + 
    Age + Gender + Income + City, family = binomial, data = negative)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0698  -0.7018  -0.6246  -0.5315   2.1232  

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)   
(Intercept)        -1.696660   0.523420  -3.241  0.00119 **
TimeYardYear       -0.001198   0.091319  -0.013  0.98953   
TimeGSYear         -0.293186   0.096818  -3.028  0.00246 **
TimeYardWeek        0.029629   0.049917   0.594  0.55281   
X.10minYardWeek     0.004619   0.059877   0.077  0.93850   
FreqVisitGS         0.015702   0.037776   0.416  0.67766   
AvgNR               0.229938   0.131053   1.755  0.07934 . 
HoursLastWeekGroup  0.008057   0.023582   0.342  0.73262   
Age                -0.043269   0.021898  -1.976  0.04816 * 
Gender              0.002840   0.147026   0.019  0.98459   
Income             -0.065472   0.026146  -2.504  0.01228 * 
CitySydney          0.129282   0.143189   0.903  0.36659   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1271.1  on 1251  degrees of freedom
Residual deviance: 1246.9  on 1240  degrees of freedom
  (420 observations deleted due to missingness)
AIC: 1270.9

Number of Fisher Scoring iterations: 4

PWI<-lm(PWI~X.10minYardWeek+TimeYardWeek+FreqVisitGS+HoursLastWeekGroup+AvgNR+
          Age+Gender+Income+City,data=data)
summary(PWI)

Call:
lm(formula = PWI ~ X.10minYardWeek + TimeYardWeek + FreqVisitGS + 
    HoursLastWeekGroup + AvgNR + Age + Gender + Income + City, 
    data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-43.994  -6.136   1.281   7.253  27.425 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        34.87270    1.83496  19.005  < 2e-16 ***
X.10minYardWeek     0.70754    0.18705   3.783 0.000160 ***
TimeYardWeek       -0.22548    0.18493  -1.219 0.222913    
FreqVisitGS         0.48300    0.12707   3.801 0.000149 ***
HoursLastWeekGroup  0.05965    0.08168   0.730 0.465299    
AvgNR               1.19358    0.45116   2.646 0.008227 ** 
Age                 0.48498    0.07618   6.366 2.47e-10 ***
Gender              0.02198    0.50953   0.043 0.965593    
Income              0.68450    0.09285   7.372 2.57e-13 ***
CitySydney         -0.86404    0.49639  -1.741 0.081916 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.37 on 1775 degrees of freedom
  (299 observations deleted due to missingness)
Multiple R-squared:  0.1056,    Adjusted R-squared:  0.101 
F-statistic: 23.27 on 9 and 1775 DF,  p-value: < 2.2e-16

library(ggplot2)
library(ggpubr)

data$FreqVisitGSF<-as.factor(data$FreqVisitGS)
data$X.10minYardWeekF<-as.factor(data$X.10minYardWeek)
a<-ggplot(data, aes(x = X.10minYardWeekF, y = PWI)) + 
  geom_boxplot(color="#69b3a2", fill="69b3a2", alpha=0.2) +
  geom_jitter(width = 0.08,col="grey",alpha=0.8)+
  labs(x = "Frequency of yard visits",y="Personal Wellbeing index")+
  theme_bw()+
  theme(axis.text.x = element_text(color="black"),
          axis.text.y = element_text(color="black"))

b<-ggplot(data, aes(x = FreqVisitGSF, y = PWI)) + 
  geom_boxplot(color="#267532", fill="#267532", alpha=0.2) +
  geom_jitter(width = 0.08,col="grey",alpha=0.8)+
  labs(x = "Frequency of public greenspace visits",y="Personal Wellbeing index")+
  theme_bw()+
  theme(axis.text.x = element_text(color="black"),
          axis.text.y = element_text(color="black"))

ggarrange(b, a, labels = c("a", "b"),
                    ncol = 2, nrow = 1)

NA
NA

a<-ggplot(positive, aes(TimeGSYear, HealthChangeYearF)) + 
    geom_jitter(position = position_jitter(width = 0.25, height = 0.02),col="grey",alpha=0.8) + 
    geom_smooth(color="#267532",method = "glm", se = T, 
        method.args = list(family = "binomial"))+
  labs(x = "Change in public greenspace visits",y="Likelihood of positive change")+
  theme_bw()+
  theme(axis.text.x = element_text(color="black"),
          axis.text.y = element_text(color="black"))

b<-ggplot(positive, aes(TimeYardYear, HealthChangeYearF)) + 
    geom_jitter(position = position_jitter(width = 0.25, height = 0.02),col="grey",alpha=0.8) + 
    geom_smooth(color="#69b3a2",method = "glm", se = T, 
        method.args = list(family = "binomial"))+
  labs(x = "Change in yard visits",y="Likelihood of positive change")+
  theme_bw()+
  theme(axis.text.x = element_text(color="black"),
          axis.text.y = element_text(color="black"))

c<-ggplot(negative, aes(TimeGSYear, HealthChangeYearF*-1)) + 
    geom_jitter(position = position_jitter(width = 0.25, height = 0.02),col="grey",alpha=0.8) + 
    geom_smooth(color="#267532",method = "glm", se = T, 
        method.args = list(family = "binomial"))+
  labs(x = "Change in public greenspace visits",y="Likelihood of negative change")+
  theme_bw()+
  theme(axis.text.x = element_text(color="black"),
          axis.text.y = element_text(color="black"))

d<-ggplot(negative, aes(TimeYardYear, HealthChangeYearF*-1)) + 
    geom_jitter(position = position_jitter(width = 0.25, height = 0.02),col="grey",alpha=0.8) + 
    geom_smooth(color="#69b3a2",method = "glm", se = T, 
        method.args = list(family = "binomial"))+
  labs(x = "Change in yard visits",y="Likelihood of negative change")+
  theme_bw()+
  theme(axis.text.x = element_text(color="black"),
          axis.text.y = element_text(color="black"))

ggarrange(a,b,c,d, labels = c("a", "b","c", "d"),
                    ncol = 2, nrow = 2)

`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'Warning: Removed 194 rows containing non-finite values (`stat_smooth()`).Warning: Removed 194 rows containing missing values (`geom_point()`).`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'Warning: Removed 200 rows containing non-finite values (`stat_smooth()`).Warning: Removed 200 rows containing missing values (`geom_point()`).