##############

#get Lewandosky dataset 
get.data=function(dset="lew") {
 if(dset=="lew") {
	#loc="http://www.bishop-hill.net/storage/LskyetalPsychSciClimate.xls"
	#download.file(loc,"temp")#,mode="wb")
	#bak=read.xls("temp",sheet=1)
	loc="http://www.climateaudit.info/data/psychology/LskyetalPsychSciClimate.csv"
	download.file(loc,"temp")
	bak=read.csv("temp")
	dim(bak) #1145 32
	return(bak)
 } #lew
 if(dset=="watts") {
	loc="d:/climate/data/psychology/Attitudes Towards Science_ Survey.csv"
	q=read.csv(loc,nrow=1,header=FALSE)
	work=read.csv(loc,skip=1,
	   colClasses=c(rep("character",3),rep("numeric",35),rep("character",11) )   )
	apply(work[,4:38],2,range)
	work$ConsensSmoke=parsew(work$ConsensSmoke)
	work$ConsensHIV=parsew(work$ConsensHIV)
	#work$InCountry=parsew(work$InCountry)
	#work$IncNeigh=parsew(bak$IncNeigh)
	work$Age=parsew(work$Age)
	work=work[,c(4:9,25:28,29:30,11:24,36:38,39:41,10,44,49)]
	work[,c(1:29,33)]=6- work[,c(1:29,33)]
	#temp= work[,c(1:29,33)]==3
	#work[,c(1:29,33)][temp]=NA
	return(work)
 }#watts
}




##used to read watts
parsew=function(x) {
	#x=bak$ConsensSmoke
	x=gsub("%","",x)
	x=gsub("?","",x)
	x=gsub("+","",x)
	x=gsub("over","",x)
	x=gsub("OVER","",x)
	x=as.numeric(x)
	reject=x>100&!is.na(x);sum(reject)
	x[reject]=NA
	#x[is.na(x)]=3
	return(x)
}


##define skeptic
define_skeptic=function(bak,dset="lew") {
  	if(dset=="lew") Break=0 else Break =1
	K=nrow(bak)
	bak$class="inconsistent"
	bak$class[bak$CO2HasNegChange>=(3+Break)& bak$CO2TempUp>=(3+Break)]="warmist" 
	bak$class[bak$CO2HasNegChange<=2& bak$CO2TempUp<=2]="skydragon"
	bak$class[bak$CO2HasNegChange<=2& bak$CO2TempUp>=3]="skeptic" #include unknown for watts
	bak$class=factor(bak$class)
	return(bak)
}


##return conspiracy counts
conspiracy.count=function(bak,ind=index,dset="lew") {
    if(dset=="lew") Break=0 ; if(dset=="watts") Break=1
    #count conspiracies
	Stat=array(NA,dim=c(17,5)) 
	row.names(Stat)=ind;
	Stat=data.frame(Stat)
	names(Stat)=c("count","inconsistent","skeptic","skydragon","warmist")
	for(j in 1:13 ) {
	 i=ind[j]
	 temp= bak[,i]>=(3+Break)
	Stat[j,]=c(sum(temp),table ( bak$class[temp])  ) 
	}
	for(j in 14:15 ) {
	 i=ind[j]
	 temp= bak[,i]<=2
	Stat[j,]=c(sum(temp),table (bak$class[temp])) 
	}
	for(j in 16:17 ) {
	 i=ind[j]
	 temp= bak[,i]<50 
	Stat[j,]=c(sum(temp,na.rm=T),table (bak$class[temp]))
	}
	return(Stat)
}


get_latent=function(bak,FUN=mean) {
	K=nrow(bak)
	U=data.frame(array(NA,dim=c(K,5)))
	names(U)=c("co2","sci","env","econ","cy")
	U$econ= apply(bak[,c(1,3,4,5,6)],1,FUN)
	U$co2 =apply(bak[,c(7:10,29)],1,FUN)
	#U$cy=apply(bak[,c(13:15,17:24,26)],1, function(x) weighted.mean(x, c(rep(.742,5),.891,.742,.742,.891, rep(.742,3 ) ) ) )
	U$cy=apply(bak[,c(13:15,17:24,26)],1, FUN)	
	g= function(y, target=c(0,25,90,99,100)   ) {
	  M=length(target)
	  v=rep(M,K)
	  for(j in (M-1):1) { 
	   v[y<target[j]]=j}
	  v=v-1
	 return(v)
	}
	Env=bak[,27:32]
	for(i in 4:6) Env[,i]=g(Env[,i])	
	U$sci=apply(bak[,c(27,28)],1,FUN)
	U$env=apply(bak[,11:12],1,FUN)
	return(U)
}




##make two classes
twoclass=function(bak,dset="lew"){ 
   if(dset=="watts") Break=1 else Break=0
	temp=bak[,1:29]> (2+Break)
	bak[,1:29][temp]=1
	bak[1:29][!temp]=0
	temp=bak[,30:32]>=80
	bak[,30:32][temp]=1
	bak[,30:32][!temp]=0
	return(bak)
}


#expand range a little
prettyf=function(x) {
	  xlim0= range(x); 
	   xlim0= c(xlim0[1]-.05*diff(xlim0), xlim0[2]+.05*diff(xlim0) )
	return(xlim0)
}
 

tableplot=function(A,r=NULL,overplot=FALSE,xlab=NULL,ylab=NULL) {
	par(bg="grey10",col="yellow",col.axis="yellow",col.lab="yellow")
	  K=nrow(A); M=ncol(A)
	  if(is.null(xlab)) xlab=names(attributes(A)$dimnames)[1]
	  if(is.null(ylab)) ylab=names(attributes(A)$dimnames)[2]
	  x= as.numeric(dimnames(A)[[1]])
	  y=as.numeric(dimnames(A)[[2]])
	  X= data.frame(x= rep(x,M), y=  c(t(array(rep(y,K),dim=c(M,K)))),count=c(A))	 
 	  xlim0=prettyf(x); ylim0=prettyf(y)
	  if(is.null(r)) r= .01
	 plot(y~x,X,pch=19,col="yellow",xlim=xlim0,ylim=ylim0,cex=r*X$count,ylab=ylab,xlab=xlab)
	 if(overplot) points(y~x,X,cex=1.5)
}




olsplot=function(Data=W,xn,yn,col="cyan",lwd=2) {
	ols=lm(Data[,yn]~Data[,xn])
	b=ols$coef; 
	a= range(Data[,xn])
	lines(a,b[1]+b[2]*a,lwd=lwd,col=col,lty=3)
}

rlmplot=function(Data=W,xn,yn,col="yellow",lwd=2) {
	fm=rlm(Data[,yn]~Data[,xn])
	b=fm$coef; a= range(Data[,xn])
	lines(a,b[1]+b[2]*a,lwd=lwd,col=col,lty=3)
}