##SCRIPT FOR BLOG POST 2012-09-08 ON LEWANDOWSKY

library(xlsReadWrite)
loc="http://www.bishop-hill.net/storage/LskyetalPsychSciClimate.xls"
	download.file(loc,"temp",mode="wb")
	bak=read.xls("temp",sheet=1)
	dim(bak) #1145 32

######
index=c(13:24,26:28,30,31)
names(bak)[index]
# [1] "CYNewWorldOrder" "CYSARS"          "CYPearlHarbor"   "CYAIDS"          "CYMLK"          
# [6] "CYMoon"          "CYArea51"        "CYJFK"           "CY911"           "CYRoswell"      
#[11] "CYDiana"         "CYOkla"          "CYCoke"          "CauseHIV"        "CauseSmoke"     
#[16] "ConsensHIV"      "ConsensSmoke"   

crankindex =c(13,14, 16, 18, 19, 21, 22, 23)
crank=apply(bak[,crankindex]>=3,1,sum)
name0=gsub("CY","",names(bak)[index])


######
## Simple Counts
########
#skydragons and skeptics
rbind( NegChange=with(bak,table(CO2HasNegChange) ),
 	TempUp=with(bak,table(CO2TempUp) ) )
#               1   2   3   4
# NegChange    133 116 356 540
# TempUp        58  69 198 820
 
with(bak,table(CO2HasNegChange,CO2TempUp) )
#              CO2TempUp
#CO2HasNegChange   1   2   3   4
#              1  53  41  36   3
#              2   3  25  43  45
#              3   0   3  99 254
#              4   2   0  20 518


#define skeptic
	temp1=bak$CO2HasNegChange<=2& bak$CO2TempUp>=3
	benchSkep=sum(temp1)/1145 #,111

#define skydragon
	temp2=bak$CO2HasNegChange<=2& bak$CO2TempUp<=2
	(benchSky=sum(temp2)/1145) # 0.1065

#define combo
	temp0=bak$CO2HasNegChange<=2
	benchComb=sum(temp0)/1145 #0.217


################
## COUNT CONSPIRACY ADHERENTS BY SUBGROUP
############################
count= count_old=c(1145,sum(temp2),sum(temp1),sum(temp0),sum(!temp1&!temp2) ) 

Stat=array(NA,dim=c(17,5)) 
row.names(Stat)=name0;
Stat=data.frame(Stat)
names(Stat)=c("count","skydragon","skeptic","combo","warmist")
for(j in 1:13 ) {
	i=index[j]
	temp= bak[,i]>=3
	Stat[j,]= c(sum(temp),sum(temp&temp2),sum(temp&temp1),sum(temp&temp0),sum(temp&!temp1&!temp2) ) 
}
for(j in 14:15 ) {
	i=index[j]
	temp= bak[,i]<3
	Stat[j,]= c(sum(temp),sum(temp&temp2),sum(temp&temp1),sum(temp&temp0),sum(temp&!temp1&!temp2) ) 
}
for(j in 16:17 ) {
	i=index[j]
	temp= bak[,i]<50 
	Stat[j,]= c(sum(temp),sum(temp&temp2),sum(temp&temp1),sum(temp&temp0),sum(temp&!temp1&!temp2) ) 
}

##############
##RESULTS
##################

	A=Stat
	row.names(A)=paste("#",row.names(A))
	A

#                count skydragon skeptic combo warmist
# NewWorldOrder    70        44       7    51      19
# SARS             42        13       3    16      26
# PearlHarbor     146        25      16    41     105
# AIDS              9         4       0     4       5
# MLK              90        13       5    18      72
# Moon             10         4       0     4       6
# Area51           35        11       2    13      22
# JFK             247        45      29    74     173
# 911              69        12       3    15      54
# Roswell          47        14       6    20      27
# Diana            25         9       1    10      15
# Okla            289        27      24    51     238
# Coke            151        24      20    44     107
# CauseHIV         16        12       1    13       3
# CauseSmoke       11         7       0     7       4
# ConsensHIV       14         8       1     9       5
# ConsensSmoke     11         4       0     4       7


############
##SCATTER PLOT 1 - COMBO
###########################

#setwd("d:/climate/data/psychology")
#png(file="lewandowsky_scatter1.png",w=600,h=600)
par(mar=c(4,4,2,1))
	Stat$share=Stat$combo/Stat$count
	plot(share~count,Stat,pch=19,xlab="Conspiracy Count",ylab="'Skeptic' Share" )
	title("'Skeptic' Share versus Conspiracy Count")

	for(i in c(1:5,7:11,13,14,15,16) ) text(Stat[i,1],Stat[i,"share"],name0[i],col=2,pos=4,font=2)
	i=6;text(Stat[i,1],Stat[i,"share"]+.01,font=2,name0[i],col=2,pos=1)
	i=12;text(Stat[i,1],Stat[i,"share"],name0[i],font=2,col=2,pos=2)
	i=17;text(Stat[i,1]+20,Stat[i,"share"],name0[i],font=2,col=2,pos=1)
	abline(h=benchComb)
#dev.off()

#####################
##SCATTER PLOT 2- LUKEWARMER
###################

#png(file="lewandowsky_scatter2.png",w=600,h=600)
par(mar=c(4,4,2,1))
	Stat$share=Stat$skeptic/Stat$count
	plot(share~count,Stat,pch=19,xlab="Conspiracy Count",ylab="'Skeptic' Share",ylim=c(0,.8) ,yaxs="i")
	title("Scatter Plot - 'True' Skeptic")

	 abline(h=benchSkep)
	#for(i in 1:17 ) text(Stat[i,1],Stat[i,"share"],name0[i],col=2,pos=4,font=2)
#dev.off()


####################
##SCATTER PLOT 3 - SKYDRAGON
#########################

#png(file="lewandowsky_scatter3.png",w=600,h=600)
par(mar=c(4,4,2,1))
	Stat$share=Stat$skydragon/Stat$count
	plot(share~count,Stat,pch=19,xlab="Conspiracy Count",ylab="'Skeptic' Share",ylim=c(0,.8) ,yaxs="i")
	title("Scatter Plot - Skydragon")
	abline(h=benchSky)
	for(i in 1:17 ) text(Stat[i,1],Stat[i,"share"],name0[i],col=2,pos=4,font=2)
#dev.off()


###########
##BARPLOT
######################
order1=order(Stat$count)

png(file="lewandowsky_barplot1.png",w=600,h=600)
par(mar=c(11,4,2,1))
  	barplot(t(Stat[order1,2:3]),beside=TRUE,las=3,ylab="Conspiracy Count",ylim=c(0,50))
 	legend("topleft",fill=c("black","lightblue"),legend=c("Skydragon","Lukewarmer") )
	box()

dev.off()


png(file="lewandowsky_barplot2.png",w=600,h=600)

work= round(t( t(Stat)/count) ,3)
 barplot( t(work[order1,c(2,3,5)]),beside=TRUE,las=3,col=1:3,ylim=c(0,.45),ylab="Fraction Adherence")
box()
	legend("topleft",fill=1:3,legend=c("Skydragon","Lukewarmer","Warmist") )
	title("Conspiracy Adherence")
dev.off()