#Identifying pertinent variables for nonresponse survey: ESS 2010

#Load the necessary packages

library(foreign)
library(descr)
library(fmsb)
library ("Hmisc")
library(hexbin)
library("psych")
library("fmsb")
library("lawstat")
library(MASS)  
library("car")
library("gmodels")
library(survey)
#Read the correct file.

setwd("I:/Work Caro/Documents/nonresponse surveys/ESS")
sink(file="es5_corrmean.txt",split=T)
library(foreign)
library(descr)
library(fmsb)
es5NR<-read.spss("ess10 for R1.sav",use.missings=T,use.value.labels=T)
es5<-as.data.frame(es5NR)
table(es5$mainsurv)

#Data cleaning and variables recoding. Recoded variable finishing with M originate from main survey
#Variable from non-response survey.

es5$satdem<-es5$q1
es5$satdem[es5$q1==11]<-NA

es5$immig<-es5$q2
es5$immig[es5$q2==11]<-NA

es5$trujus<-es5$q3
es5$trujus[es5$q3==11]<-NA

es5$accpol<-es5$q4
es5$accpol[es5$q4==11]<-NA

es5$polintr<-5-es5$q5
es5$polintr[es5$q5==5]<-NA

es5$survey<-6-es5$q6_a
es5$survey[es5$q6_a==6]<-NA

es5$science<-6-es5$q6_b
es5$science[es5$q6_b==6]<-NA

es5$happy<-es5$q7
es5$happy[es5$q7==11]<-NA

es5$health<-6-es5$q8
es5$health[es5$q8==6]<-NA

es5$socmeet<-es5$q9
es5$socmeet[es5$q9==8]<-NA


es5$socac<-es5$q10
es5$socac[es5$q10==6]<-NA


es5$safe<-5-es5$q11
es5$safe[es5$q11==5]<-NA


es5$educr<-es5$q12
es5$educr[es5$q12<1|es5$q12>9]<-NA


es5$work[es5$q13>=2&es5$q13<7]<-0
es5$work[es5$q13==1]<-1
es5$work[es5$q13==7]<-NA

es5$tired<-NA
es5$tired<-es5$q14_a
es5$tired[es5$q14_a==6]<-NA


es5$worry<-NA
es5$worry<-es5$q14_b
es5$worry[es5$q14_b==6]<-NA
table(es5$worry)

es5$fix[es5$q15==2]<-0
es5$fix[es5$q15==1]<-1
es5$fix[es5$q15>=3]<-NA


es5$mobile<-NA
es5$mobile[es5$q16==2]<-0
es5$mobile[es5$q16==1]<-1

es5$dir<-NA
es5$dir[es5$q15b==2 & es5$q16b==2]<-0
es5$dir[es5$q15b==1|es5$q16b==1]<-1
es5$dir[es5$q15b>=3&es5$q16b>=3]<-NA

es5$alone<-NA
es5$alone[es5$q17==1]<-1
es5$alone[es5$q17>1]<-0
es5$alone[es5$q17>=90]<-NA


es5$gndr<-NA
es5$gndr[es5$q19==2]<-0
es5$gndr[es5$q19==1]<-1
es5$gndr[es5$q19>=3]<-NA

es5$ager<-2010-es5$q19b
es5$ager[es5$q19b==9999]<-NA


#variable from main survey.

es5$satdemM<-es5$STFDEM
es5$satdemM[es5$STFDEM>=77]<-NA


es5$immigM<-es5$IMWBCNT
es5$immigM[es5$IMWBCNT>=77]<-NA

es5$trujusM<-es5$TRSTLGL
es5$trujusM[es5$TRSTLGL>=77]<-NA

es5$accpolM<-es5$BPLCDC
es5$accpolM[es5$BPLCDC>=77]<-NA

es5$polintrM<-es5$POLINTR
es5$polintrM[es5$POLINTR>=7]<-NA

es5$surveyM<-es5$S1qch
es5$surveyM[es5$S1qch==8]<-NA

es5$scienceM<-es5$SCNSENV
es5$scienceM[es5$SCNSENV>=7]<-NA

es5$happyM<-es5$HAPPY
es5$happyM[es5$HAPPY>=77]<-NA


es5$healthM<-es5$HEALTH
es5$healthM[es5$HEALTH>=7]<-NA


es5$socmeetM<-es5$SCLMEET
es5$socmeetM[es5$SCLMEET>=8]<-NA

es5$socacM<-es5$SCLACT
es5$socacM[es5$SCLACT>=7]<-NA

es5$safeM<-es5$AESFDRK
es5$safeM[es5$AESFDRK>=7]<-NA

es5$educrM<-es5$EDLVDCH_c
es5$educrM[es5$EDLVDCH_c<1|es5$EDLVDCH_c>9]<-1
es5$educrM[es5$mainsurv==0]<-NA


es5$workM<-NA
es5$workM[es5$MAINACT_c>1]<-0
es5$workM[es5$MAINACT_c==1]<-1
es5$workM[es5$MAINACT_c>=7]<-NA

es5$tiredM<-NA
es5$tiredM<-es5$TRDAWRK
es5$tiredM[es5$TRDAWRK>=6]<-NA

es5$worryM<-es5$WRYWPRB
es5$worryM[es5$WRYWPRB>=6]<-NA

es5$fixM<-NA
es5$fixM[es5$FXLTPH==2]<-0
es5$fixM[es5$FXLTPH==1]<-1
es5$fixM[es5$FXLTPH>=7]<-NA

es5$mobileM<-NA
es5$mobileM[es5$Q4qch==2]<-0
es5$mobileM[es5$Q4qch==1]<-1
es5$mobileM[es5$Q4qch>=7]<-NA

es5$dirM<-NA
es5$dirM[es5$Q3qch==2]<-0
es5$dirM[es5$Q3qch==1]<-1
es5$dirM[es5$Q3qch>=3]<-NA

es5$aloneM<-NA
es5$aloneM[es5$HHMMB==1]<-1
es5$aloneM[es5$HHMMB>1]<-0
es5$aloneM[es5$HHMMB>=80]<-NA

es5$gndrN[es5$GNDR==2]<-0
es5$gndrN[es5$GNDR==1]<-1
es5$gndrN[es5$GNDR>=3]<-NA

es5$agerM<-2010-es5$YRBRN
es5$agerM[es5$YRBRN==9999]<-NA



#select the cases for which the respondent answered both survey and sex and year of birth corresponds.

nrandmain <- which(es5$nrepsurv==1&es5$mainsurv==1)
es51<-es5[nrandmain,]
table(es51$GNDR)
table(es51$q19)
es51$test<-es51$ager-es51$agerM
es51$agetest<-0
es51$agetest[es51$test<=-3|es51$test>=3]<-1
table(es51$agetest)
nrandmain1<-which((es51$GNDR==es51$q19| es51$q19 == 3)&es51$agetest==0)
es5corr<-es51[nrandmain1,]


type<- c("corr","mean")
es5corrmean<-as.data.frame(type)


#create difference varaibles for mean.correlation 
es5corr$satdemD<-NA
es5corr$satdemD<-es5corr$satdem-es5corr$satdemM
es5corrmean$satdem[es5corrmean$type=="corr"]<-cor(es5corr$satdem,es5corr$satdemM,"pairwise.complete.obs", method="spearman")
es5corrmean$satdem[es5corrmean$type=="mean"]<-mean(es5corr$satdemD, na.rm=TRUE)

es5corr$immigD<-NA
es5corr$immigD<-es5corr$immig-es5corr$immigM
es5corrmean$immig[es5corrmean$type=="corr"]<-cor(es5corr$immig,es5corr$immigM,"pairwise.complete.obs", method="spearman")
es5corrmean$immig[es5corrmean$type=="mean"]<-mean(es5corr$immigD, na.rm=TRUE)

es5corr$trujusD<-NA
es5corr$trujusD<-es5corr$trujus-es5corr$trujusM
es5corrmean$trujus[es5corrmean$type=="corr"]<-cor(es5corr$trujus,es5corr$trujusM,"pairwise.complete.obs", method="spearman")
es5corrmean$trujus[es5corrmean$type=="mean"]<-mean(es5corr$trujusD, na.rm=TRUE)

es5corr$accpolD<-NA
es5corr$accpolD<-es5corr$accpol-es5corr$accpolM
es5corrmean$accpol[es5corrmean$type=="corr"]<-cor(es5corr$accpol,es5corr$accpolM,"pairwise.complete.obs", method="spearman")
es5corrmean$accpol[es5corrmean$type=="mean"]<-mean(es5corr$accpolD, na.rm=TRUE)


es5corr$polintrD<-NA
es5corr$polintrD<-es5corr$polintr-es5corr$polintrM
es5corrmean$polintr[es5corrmean$type=="corr"]<-cor(es5corr$polintr,es5corr$polintrM,"pairwise.complete.obs", method="spearman")
es5corrmean$polintr[es5corrmean$type=="mean"]<-mean(es5corr$polintrD, na.rm=TRUE)

es5corr$surveyD<-NA
es5corr$surveyD<-es5corr$survey-es5corr$surveyM
es5corrmean$survey[es5corrmean$type=="corr"]<-cor(es5corr$survey,es5corr$surveyM,"pairwise.complete.obs", method="spearman")
es5corrmean$survey[es5corrmean$type=="mean"]<-mean(es5corr$surveyD, na.rm=TRUE)

es5corr$scienceD<-NA
es5corr$scienceD<-es5corr$science-es5corr$scienceM
es5corrmean$science[es5corrmean$type=="corr"]<-cor(es5corr$science,es5corr$scienceM,"pairwise.complete.obs", method="spearman")
es5corrmean$science[es5corrmean$type=="mean"]<-mean(es5corr$scienceD, na.rm=TRUE)

es5corr$happyD<-NA
es5corr$happyD<-es5corr$happy-es5corr$happyM
es5corrmean$happy[es5corrmean$type=="corr"]<-cor(es5corr$happy,es5corr$happyM,"pairwise.complete.obs", method="spearman")
es5corrmean$happy[es5corrmean$type=="mean"]<-mean(es5corr$happyD, na.rm=TRUE)
nrandmain <- which(es5corr$happyM>=0&es5corr$happy>=0)

es5corr$healthD<-NA
es5corr$healthD<-es5corr$health-es5corr$healthM
es5corrmean$health[es5corrmean$type=="corr"]<-cor(es5corr$health,es5corr$healthM,"pairwise.complete.obs", method="spearman")
es5corrmean$health[es5corrmean$type=="mean"]<-mean(es5corr$healthD, na.rm=TRUE)

es5corr$socacD<-NA
es5corr$socacD<-es5corr$socac-es5corr$socacM
es5corrmean$socac[es5corrmean$type=="corr"]<-cor(es5corr$socac,es5corr$socacM,"pairwise.complete.obs", method="spearman")
es5corrmean$socac[es5corrmean$type=="mean"]<-mean(es5corr$socacD, na.rm=TRUE)
nrandmain <- which(es5corr$socacM>=0&es5corr$socac>=0)
es5temp<-es5corr[nrandmain,]  
es5corrmean$socac[es5corrmean$type=="Cronbach alpha"]<-CronbachAlpha(data.frame(es5temp$socac,es5temp$socacM))

es5corr$socmeetD<-NA
es5corr$socmeetD<-es5corr$socmeet-es5corr$socmeetM
es5corrmean$socmeet[es5corrmean$type=="corr"]<-cor(es5corr$socmeet,es5corr$socmeetM,"pairwise.complete.obs", method="spearman")
es5corrmean$socmeet[es5corrmean$type=="mean"]<-mean(es5corr$socmeetD, na.rm=TRUE)

es5corr$safeD<-NA
es5corr$safeD<-es5corr$safe-es5corr$safeM
es5corrmean$safe[es5corrmean$type=="corr"]<-cor(es5corr$safe,es5corr$safeM,"pairwise.complete.obs", method="spearman")
es5corrmean$safe[es5corrmean$type=="mean"]<-mean(es5corr$safeD, na.rm=TRUE)

es5corr$educrD<-NA
es5corr$educrD<-es5corr$educr-es5corr$educrM
es5corrmean$educr[es5corrmean$type=="corr"]<-cor(es5corr$educr,es5corr$educrM,"pairwise.complete.obs", method="spearman")
es5corrmean$educr[es5corrmean$type=="mean"]<-mean(es5corr$educrD, na.rm=TRUE)


es5corr$workD<-NA
es5corr$workD<-es5corr$work-es5corr$workM
es5corrmean$work[es5corrmean$type=="corr"]<-phi(table(es5corr$work,es5corr$workM))
es5corrmean$work[es5corrmean$type=="mean"]<-mean(es5corr$workD, na.rm=TRUE)

es5corr$tiredD<-NA
es5corr$tiredD<-es5corr$tired-es5corr$tiredM
table(es5corr$tired,es5corr$tiredM)
es5corrmean$tired[es5corrmean$type=="corr"]<-cor(es5corr$tired,es5corr$tiredM,"pairwise.complete.obs", method="spearman")
es5corrmean$tired[es5corrmean$type=="mean"]<-mean(es5corr$tiredD, na.rm=TRUE)

es5corr$worryD<-NA
es5corr$worryD<-es5corr$worry-es5corr$worryM
es5corrmean$worry[es5corrmean$type=="corr"]<-cor(es5corr$worry,es5corr$worryM,"pairwise.complete.obs", method="spearman")
es5corrmean$worry[es5corrmean$type=="mean"]<-mean(es5corr$worryD, na.rm=TRUE)

es5corr$fixD<-NA
es5corr$fixD<-es5corr$fix-es5corr$fixM
es5corrmean$fix[es5corrmean$type=="corr"]<-phi(table(es5corr$fix,es5corr$fixM))
es5corrmean$fix[es5corrmean$type=="mean"]<-mean(es5corr$fixD, na.rm=TRUE)


es5corr$dirD<-NA
es5corr$dirD<-es5corr$dir-es5corr$dirM
es5corrmean$dir[es5corrmean$type=="corr"]<-phi(table(es5corr$dir,es5corr$dirM))
es5corrmean$dir[es5corrmean$type=="mean"]<-mean(es5corr$dirD, na.rm=TRUE)


es5corr$mobileD<-NA
es5corr$mobileD<-es5corr$mobile-es5corr$mobileM
es5corrmean$mobile[es5corrmean$type=="corr"]<-phi(table(es5corr$mobile,es5corr$mobileM))
es5corrmean$mobile[es5corrmean$type=="mean"]<-mean(es5corr$mobileD, na.rm=TRUE)

es5corr$gndrD<-NA
es5corr$gndrD<-es5corr$gndr-es5corr$gndrN
es5corrmean$gndr[es5corrmean$type=="corr"]<-phi(table(es5corr$gndr,es5corr$gndrN))
es5corrmean$gndr[es5corrmean$type=="mean"]<-mean(es5corr$gndrD, na.rm=TRUE)

es5corr$agerD<-NA
es5corr$agerD<-es5corr$ager-es5corr$agerM
es5corrmean$ager[es5corrmean$type=="corr"]<-cor(es5corr$ager,es5corr$agerM,"pairwise.complete.obs", method="spearman")
es5corrmean$ager[es5corrmean$type=="mean"]<-mean(es5corr$agerD, na.rm=TRUE)
table(es5corr$agerD)

es5corr$aloneD<-NA
es5corr$aloneD<-es5corr$alone-es5corr$aloneM
es5corrmean$alone[es5corrmean$type=="corr"]<-phi(table(es5corr$alone,es5corr$aloneM))
es5corrmean$alone[es5corrmean$type=="mean"]<-mean(es5corr$aloneD, na.rm=TRUE)

boxplot(es5corr$satdemD,es5corr$polintrD, es5corr$immigD,es5corr$trujusD,es5corr$accpolD,es5corr$surveyD,es5corr$scienceD,es5corr$happyD,es5corr$healthD,es5corr$socacD,
        es5corr$socmeetD,es5corr$safeD,es5corr$educrD,es5corr$tiredD,es5corr$worryD,es5corr$ageD, names=c("de","pi","im","tj","plc","su","sc","ha","he","sa","sm"
                                                                                                          ,"sf","ed","ti","wo","ag"))
boxplot(es5corr$satdemD,es5corr$polintrD, es5corr$immigD,es5corr$trujusD,es5corr$healthD,es5corr$socacD,es5corr$socmeetD,es5corr$agerD,
        names=c("demo","polin","immig","trstj","heal","soca","socm","age"))  

#t-test and Wilcoxon test for differencec in distributions
 

wilcox.test(es5corr$satdem,es5corr$satdemM,paired=TRUE)
wilcox.test(es5corr$polintr,es5corr$polintrM,paired=TRUE)
wilcox.test(es5corr$immig,es5corr$immigM,paired=TRUE)
wilcox.test(es5corr$trujus,es5corr$trujusM,paired=TRUE)
wilcox.test(es5corr$accpol,es5corr$accpolM,paired=TRUE)
wilcox.test(es5corr$survey,es5corr$surveyM,paired=TRUE)
wilcox.test(es5corr$science,es5corr$scienceM,paired=TRUE)
wilcox.test(es5corr$happy,es5corr$happyM,paired=TRUE)
wilcox.test(es5corr$healthM,es5corr$health,paired=TRUE)
wilcox.test(es5corr$socac,es5corr$socacM,paired=TRUE)
wilcox.test(es5corr$socmeet,es5corr$socmeetM,paired=TRUE)
wilcox.test(es5corr$safe,es5corr$safeM,paired=TRUE)
wilcox.test(es5corr$educr,es5corr$educrM,paired=TRUE)
wilcox.test(es5corr$tired,es5corr$tiredM,paired=TRUE)
wilcox.test(es5corr$worry,es5corr$worryM,paired=TRUE)
t.test(es5corr$ager,es5corr$agerM,paired=TRUE)


table(es5corr$work)
table(es5corr$workM)
prop.test(matrix(c(125,120,132,140),nrow=2,ncol=2))
table(es5corr$alone)
table(es5corr$aloneM)
prop.test(matrix(c(208,214,37,46),nrow=2,ncol=2))
table(es5corr$fix)
table(es5corr$fixM)
prop.test(matrix(c(23,17,239,249),nrow=2,ncol=2))
table(es5corr$mobile)
table(es5corr$mobileM)
prop.test(matrix(c(36,39,226,223),nrow=2,ncol=2))
table(es5corr$dir)
table(es5corr$dirM)
prop.test(matrix(c(13,10,211,235),nrow=2,ncol=2))


#creation of harmonized variable: takes the values of the main survey for respondents and of NRS for nonrespondents.
es5<-es5[order(-es5$mainsurv,-es5$nrepsurv),]
v1<-es5$agerM[es5$mainsurv==1]
v2<-es5$ager[es5$mainsurv==0]
es5$ager_h<-c(v1,v2)
table(es5$count)
es5$agecat<-NA
es5$agecat[es5$ager_h<30]<-1
es5$agecat[es5$ager_h>=30&es5$ager_h<=65]<-2
es5$agecat[es5$ager_h>65]<-3

v1<-es5$educrM[es5$mainsurv==1]
v2<-es5$educr[es5$mainsurv==0]
es5$educr_h<-c(v1,v2)
table(es5$educr_h)
es5$educat<-NA
es5$educat[es5$educr_h<=2]<-0
es5$educat[es5$educr_h>2]<-1
table(es5$educr_h, es5$mainsurv)

v1<-es5$aloneM[es5$mainsurv==1]
v2<-es5$alone[es5$mainsurv==0]
es5$alone_h<-c(v1,v2)
table(es5$alone_h)

v1<-es5$gndrN[es5$mainsurv==1]
v2<-es5$gndr[es5$mainsurv==0]
es5$gndr_h<-c(v1,v2)
table(es5$gndr_h)

v1<-es5$agerM[es5$mainsurv==1]
v2<-es5$ager[es5$mainsurv==0]
es5$ager_h<-c(v1,v2)
table(es5$ager_h)

v1<-es5$educrM[es5$mainsurv==1]
v2<-es5$educr[es5$mainsurv==0]
es5$educr_h<-c(v1,v2)
table(es5$educr_h)

v1<-es5$workM[es5$mainsurv==1]
v2<-es5$work[es5$mainsurv==0]
es5$work_h<-c(v1,v2)
table(es5$work_h)
v1<-es5$polintrM[es5$mainsurv==1]
v2<-es5$polintr[es5$mainsurv==0]
es5$polintr_h<-c(v1,v2)
table(es5$polintr_h)

v1<-es5$immigM[es5$mainsurv==1]
v2<-es5$immig[es5$mainsurv==0]
es5$immig_h<-c(v1,v2)

v1<-es5$surveyM[es5$mainsurv==1]
v2<-es5$survey[es5$mainsurv==0]
es5$survey_h<-c(v1,v2)

v1<-es5$scienceM[es5$mainsurv==1]
v2<-es5$science[es5$mainsurv==0]
es5$science_h<-c(v1,v2)


v1<-es5$happyM[es5$mainsurv==1]
v2<-es5$happy[es5$mainsurv==0]
es5$happy_h<-c(v1,v2)

v1<-es5$worryM[es5$mainsurv==1]
v2<-es5$worry[es5$mainsurv==0]
es5$worried_h<-c(v1,v2)
v1<-es5$tiredM[es5$mainsurv==1]
v2<-es5$tired[es5$mainsurv==0]
es5$tired_h<-c(v1,v2)
v1<-es5$healthM[es5$mainsurv==1]
v2<-es5$health[es5$mainsurv==0]
es5$health_h<-c(v1,v2)

v1<-es5$socmeetM[es5$mainsurv==1]
v2<-es5$socmeet[es5$mainsurv==0]
es5$socmeet_h<-c(v1,v2)


v1<-es5$satdemM[es5$mainsurv==1]
v2<-es5$satdem[es5$mainsurv==0]
es5$satdem_h<-c(v1,v2)
table(es5$satdem_h)
es5$satdemcat<-NA
es5$satdemcat[es5$satdem_h>=7]<-1
es5$satdemcat[es5$satdem_h<7]<-0

v1<-es5$socacM[es5$mainsurv==1]
               
v2<-es5$socac[es5$mainsurv==0]
              
              
es5$socac_h<-c(v1,v2)
table(es5$socac_h,es5$mainsurv)

v1<-es5$fixM[es5$mainsurv==1]
v2<-es5$fix[es5$mainsurv==0]
es5$fix_h<-c(v1,v2)
table(es5$fix_h)

v1<-es5$mobileM[es5$mainsurv==1]
v2<-es5$mobile[es5$mainsurv==0]
es5$mobile_h<-c(v1,v2)
table(es5$mobile_h)

v1<-es5$dirM[es5$mainsurv==1]
v2<-es5$dir[es5$mainsurv==0]
es5$dir_h<-c(v1,v2)
table(es5$dir_h)

v1<-es5$trujusM[es5$mainsurv==1]
v2<-es5$trujus[es5$mainsurv==0]
es5$trujus_h<-c(v1,v2)
table(es5$trujus_h)
table(v1)

es5$trujuscat<-NA
es5$trujuscat[es5$trujus_h<4]<-0
es5$trujuscat[es5$trujus_h>=4&es5$trujus_h<=7]<-1
es5$trujuscat[es5$trujus_h>7]<-2


#compare the distrubutions for respondents and non-respondents.
crosstab(es5$gndr_h,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$agecat,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$educat,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$work_h,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$alone_h,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$polintr_h,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$satdemcat,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$satdem_h,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$socac_h,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$fix_h,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$mobile_h,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$dir_h,es5$mainsurv,prop.c=T,chisq=T)
crosstab(es5$trujuscat,es5$mainsurv,prop.c=T,chisq=T)




#review
es5NR<-es5[which(es5$nrepsurv==1),]
es5NR<-es5NR[which(!is.na(es5NR$ager)&!is.na(es5NR$gndr)&!is.na(es5NR$educr)&!is.na(es5NR$work)&!is.na(es5NR$alone)&!is.na(es5NR$dir)
                   &!is.na(es5NR$mobile)&!is.na(es5NR$fix)&!is.na(es5NR$socac)&!is.na(es5NR$socmeet)&!is.na(es5NR$health)&!is.na(es5NR$happy)
                   &!is.na(es5NR$science)&!is.na(es5NR$survey)&!is.na(es5NR$polintr)&!is.na(es5NR$trujus)
                   &!is.na(es5NR$satdem) &!is.na(es5NR$immig)),]
es5NR$NR<-1
es5NR$NR[es5NR$mainsurv==1]<-0
table(es5NR$NR)
es5review<-glm(NR~ager+gndr+educr+work+alone+fix+mobile+dir+socac+socmeet+health+happy
                +science+survey+polintr+trujus+immig+satdem,data=es5NR,family=binomial)
summary(es5review)
NagelkerkeR2(es5review)

es5NR$propscore<-predict(es5review, type="response")
es5NR$propweight<-1/es5NR$propscore

type<- c("RUmean","RWmean", "NRmean")
es5r<-as.data.frame(type)
es5NRNR<-es5NR[which(es5NR$NR==1),]

es5r$satdem[es5r$type=="RUmean"]<-mean(es5NR$satdemM, na.rm=TRUE)
es5r$satdem[es5r$type=="RWmean"]<-weighted.mean(es5NR$satdemM,es5NR$propweight, na.rm=TRUE)
es5r$satdem[es5r$type=="NRmean"]<-mean(es5NR$satdem, na.rm=TRUE)

es5r$immig[es5r$type=="RUmean"]<-mean(es5NR$immigM, na.rm=TRUE)
es5r$immig[es5r$type=="RWmean"]<-weighted.mean(es5NR$immigM,es5NR$propweight, na.rm=TRUE)
es5r$immig[es5r$type=="NRmean"]<-mean(es5NR$immig, na.rm=TRUE)

es5r$trujus[es5r$type=="RUmean"]<-mean(es5NR$trujusM, na.rm=TRUE)
es5r$trujus[es5r$type=="RWmean"]<-weighted.mean(es5NR$trujusM,es5NR$propweight, na.rm=TRUE)
es5r$trujus[es5r$type=="NRmean"]<-mean(es5NR$trujus, na.rm=TRUE)

es5r$polintr[es5r$type=="RUmean"]<-mean(es5NR$polintrM, na.rm=TRUE)
es5r$polintr[es5r$type=="RWmean"]<-weighted.mean(es5NR$polintrM,es5NR$propweight, na.rm=TRUE)
es5r$polintr[es5r$type=="NRmean"]<-mean(es5NR$polintr, na.rm=TRUE)

es5r$survey[es5r$type=="RUmean"]<-mean(es5NR$surveyM, na.rm=TRUE)
es5r$survey[es5r$type=="RWmean"]<-weighted.mean(es5NR$surveyM,es5NR$propweight, na.rm=TRUE)
es5r$survey[es5r$type=="NRmean"]<-mean(es5NR$survey, na.rm=TRUE)

es5r$science[es5r$type=="RUmean"]<-mean(es5NR$scienceM, na.rm=TRUE)
es5r$science[es5r$type=="RWmean"]<-weighted.mean(es5NR$scienceM,es5NR$propweight, na.rm=TRUE)
es5r$science[es5r$type=="NRmean"]<-mean(es5NR$science, na.rm=TRUE)

es5r$happy[es5r$type=="RUmean"]<-mean(es5NR$happyM, na.rm=TRUE)
es5r$happy[es5r$type=="RWmean"]<-weighted.mean(es5NR$happyM,es5NR$propweight, na.rm=TRUE)
es5r$happy[es5r$type=="NRmean"]<-mean(es5NR$happy, na.rm=TRUE)

es5r$health[es5r$type=="RUmean"]<-mean(es5NR$healthM, na.rm=TRUE)
es5r$health[es5r$type=="RWmean"]<-weighted.mean(es5NR$healthM,es5NR$propweight, na.rm=TRUE)
es5r$health[es5r$type=="NRmean"]<-mean(es5NR$health, na.rm=TRUE)


es5r$socmeet[es5r$type=="RUmean"]<-mean(es5NR$socmeetM, na.rm=TRUE)
es5r$socmeet[es5r$type=="RWmean"]<-weighted.mean(es5NR$socmeetM,es5NR$propweight, na.rm=TRUE)
es5r$socmeet[es5r$type=="NRmean"]<-mean(es5NR$socmeet, na.rm=TRUE)

es5r$socac[es5r$type=="RUmean"]<-mean(es5NR$socacM, na.rm=TRUE)
es5r$socac[es5r$type=="RWmean"]<-weighted.mean(es5NR$socacM,es5NR$propweight, na.rm=TRUE)
es5r$socac[es5r$type=="NRmean"]<-mean(es5NR$socac, na.rm=TRUE)

es5r$tired[es5r$type=="RUmean"]<-mean(es5NR$tiredM, na.rm=TRUE)
es5r$tired[es5r$type=="RWmean"]<-weighted.mean(es5NR$tiredM,es5NR$propweight, na.rm=TRUE)
es5r$tired[es5r$type=="NRmean"]<-mean(es5NR$tired, na.rm=TRUE)

es5r$worry[es5r$type=="RUmean"]<-mean(es5NR$worryM, na.rm=TRUE)
es5r$worry[es5r$type=="RWmean"]<-weighted.mean(es5NR$worryM,es5NR$propweight, na.rm=TRUE)
es5r$worry[es5r$type=="NRmean"]<-mean(es5NR$worry, na.rm=TRUE)

write.table(es5r,file= "es5r.txt",sep="\t") 

#significance
#create the right dataset
es5bothw<-es5NR[which(es5NR$mainsurv==1),]
es5bothw$propweight2<-es5bothw$propweight
es5bothw$type<-1
es5bothw$type2<-1

es5bothnw<-es5bothw
es5bothnw$propweight<-1
es5bothnw$propweight2<-0
es5bothnw$type<-0
es5bothnw$type2<-NA

es5NRnw<-es5NR
es5NRnw$propweight<-0
es5NRnw$propweight2<-1
es5NRnw$type<-NA
es5NRnw$type2<-0

es5sign<-rbind(es5bothw,es5bothnw,es5NRnw)
table(es5sign$type2)

design1<-svydesign(ids=~1,weights=es5sign$propweight, data=es5sign)

svychisq(~type+satdem_h, design = design1, statistic = "F") 
svychisq(~type+immig_h, design = design1, statistic = "F") 
svychisq(~type+trujus_h, design = design1, statistic = "F") 
svychisq(~type+polintr_h, design = design1, statistic = "F")
svychisq(~type+survey_h, design = design1, statistic = "F")
svychisq(~type+science_h, design = design1, statistic = "F")
svychisq(~type+happy_h, design = design1, statistic = "F")
svychisq(~type+health_h, design = design1, statistic = "F")
svychisq(~type+socmeet_h, design = design1, statistic = "F")
svychisq(~type+socac_h, design = design1, statistic = "F")
svychisq(~type+worried_h, design = design1, statistic = "F")
svychisq(~type+tired_h, design = design1, statistic = "F")
#all p-value above .95
design2<-svydesign(ids=~1,weights=es5sign$propweight2, data=es5sign)
svymean(es5sign$satdem_h,design=design2, na.rm=TRUE)
svychisq(~type2+immig_h, design = design2, statistic = "F") 
svychisq(~type2+trujus_h, design = design2, statistic = "F") 
svychisq(~type2+polintr_h, design = design2, statistic = "F") 
svychisq(~type2+survey_h, design = design2, statistic = "F")
svychisq(~type2+science_h, design = design2, statistic = "F")
svychisq(~type2+happy_h, design = design2, statistic = "F")
svychisq(~type2+socmeet_h, design = design2, statistic = "F")
svychisq(~type2+health_h, design = design2, statistic = "F")
svychisq(~type2+socac_h, design = design2, statistic = "F")
svychisq(~type2+worried_h, design = design2, statistic = "F")
svychisq(~type2+tired_h, design = design2, statistic = "F")