How scientists perceive the evolutionary origin of human traits: results of a survey study

This page is a study. The page identifier is Op_en7839
Moderator:Jouni (see all)

Upload data {{#opasnet_base_link:Op_en7839}}

How scientists perceive the evolutionary origin of human traits: results of a survey study is a research study that was published in journal Ecology and Evolution in 2018.

Article

Name: How scientists perceive the evolutionary origin of human traits: results of a survey study
Authors: Hanna Tuomisto, Matleena Tuomisto, Jouni T. Tuomisto
Journal: Ecology and Evolution
Publishing year: 2018 (accepted Jan 2, 2018)
DOI: doi:10.1002/ece3.3887
Appendix 2: Human evolution report survey1
Survey data on Dryad
Survey data on Opasnet

Materials and methods

Data

This code gets the files and makes them publicly available.

+ Show code - Hide code

# This is code Op_en7839/publishdata on page [[How scientists perceive the evolutionary origin of human traits: results of a survey study]]

library(OpasnetUtils)
objects.get(token)

melted_survey$Country <- NA # To protect anonymity, country information is removed
survey1$Country <- NA

objects.store(critique, melted_survey, survey1, survey1_questions, aahhypos, allhypos, envhypos, sochypos)
cat("Data.frames critique, melted_survey, survey1, and survey1_questions; and info vectors aahhypos, allhypos, envhypos, sochypos stored.\n")

Use this code in statistical software R to download the data for further use:

install.packages("OpasnetUtils") # You need this line only if you haven't installed this package yet.
library(OpasnetUtils)
objects.latest("Op_en7839", code_name="publishdata") # This line downloads data.frames with the data

Human evolution report survey1 is a report that shows all questions in the survey and simple histograms or frequencies for each question.

Calculations

Preprocess for Dryad data

The Dryad data does not contain objects melted_survey nor critique. They can be calculated from survey1 using this code. Note! This code does not fetch anything, so it only works on your own computer if you have already downloaded object survey1.

+ Show code - Hide code

# This is code Op_en7839/ on page [[How scientists perceive the evolutionary origin of human traits: results of a survey study]]

melted_survey <- melt(
  survey1[-c(5:22,80:123)], # Remove columns with raw expertise answers and about AAH criticism etc.
  measure.vars = allhypos, 
  variable.name = "Hypothesis",
  value.name="Answer"
)

melted_survey$Answer <- factor(
  melted_survey$Answer,
  levels = c(
    "Very likely",
    "Moderately likely",
    "No opinion",
    "Moderately unlikely",
    "Very unlikely"
  ),
  ordered = TRUE
)

# Calculates numbers of likely answers for each hypothesis
ord <- aggregate(melted_survey$Answer %in% c("Very likely", "Moderately likely"), melted_survey["Hypothesis"], sum) 

# Orders Hypothesis column to factor with most popular first
melted_survey$Hypothesis <- factor(melted_survey$Hypothesis, levels=ord$Hypothesis[order(-ord$x)], ordered = TRUE)

melted_survey <- merge(
  melted_survey, 
  survey1_questions[c("Name", "Trait")], # Remove Obs and Question columns, keep hypothesis and trait.
  by.x="Hypothesis",
  by.y="Name"
)

melted_survey$Trait <- factor(
  melted_survey$Trait,
  levels =  c(levels(melted_survey$Trait)[levels(melted_survey$Trait) != "Other"], "Other"),
  ordered = TRUE
)

levels(melted_survey$Hypothesis) <- gsub("\\.", " ", levels(melted_survey$Hypothesis))

################ AAH critique

critique <- melt(
  survey1[c(125, 111 ,89:108)],
  id.vars = 1:2,
  variable.name = "Critique",
  value.name = "Answer"
)

critique$Answer <- factor(critique$Answer, levels = c(
  "Fully agree",
  "Mostly agree",
  "No opinion",
  "Mostly disagree",
  "Strongly disagree",
  NA
), ordered = TRUE
)

ord <- aggregate(critique$Answer %in% c("Fully agree", "Mostly agree"), by = critique["Critique"], FUN = sum)
critique$Critique <- factor(critique$Critique, levels = ord$Critique[order(ord$x)], ordered = TRUE)

levels(critique$Critique) <- gsub("\\.", " ", levels(critique$Critique))

Correlation and histograms

Model run 11.9.2016

+ Show code - Hide code

# This is code Op_en7839/ on page [[How scientists perceive the evolutionary origin of human traits: results of a survey study]]

library(OpasnetUtils)
library(ggplot2)
library(reshape2)

objects.latest("Op_en7839", code_name="publishdata") # This line downloads data.frames with the data

#### Correlations

if(FALSE){

survey_correlations <- (cor(survey1[hypocols1], method="spearman", use="pairwise.complete.obs"))

oprint(survey_correlations)

temp <- colnames(survey_correlations)

melted_correlations <- melt(survey_correlations)

melted_correlations$Var1 <- factor(melted_correlations$Var1, levels=temp)
melted_correlations$Var2 <- factor(melted_correlations$Var2, levels=temp)
melted_correlations$value <- ifelse(melted_correlations$value >= 0.99,NA,melted_correlations$value)

ggplot(melted_correlations, aes(x = Var1, y = Var2, fill = value, label= round(value, 2)))+
  geom_raster()+
  theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4))+
  scale_fill_gradient2(low = "#480610", mid = "#FFFFFF", high = "#06480F", midpoint = 0, space = "Lab", guide = "colourbar")

}#if(FALSE)


#### Histograms

if(FALSE){

for(i in unique(melted_survey$Trait)){
  trait <- ggplot(melted_survey[melted_survey$Trait == i,],		
                  aes(Answer)) +
    geom_histogram(fill="red") +
    labs(x="", y="Count") +
    labs(title= paste(i, "in survey")) +
    theme_grey(base_size=24) +
    theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) +
    facet_wrap(~Hypothesis)
  
  print(trait)
  
}

} # if(FALSE)

survey_graph <- ggplot(melted_survey, aes(Hypothesis, fill=Answer)) +		## Plotted by trait
  geom_histogram() +
  labs(x="", y="Count", title="Everyone's views") +
  theme_grey(base_size=24) +
  theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4))

survey_graph + facet_wrap(~Trait, ncol=4, scales="free_x")

if(FALSE){

### The code to do the publishable version of the graph above.

varit <- c("mediumseagreen", "greenyellow", "gold", "darkorange", "firebrick1")
varit2 <- c("dodgerblue2", "slateblue", "darkorchid", "violetred", "brown1")

pdf(file="all_answers.pdf", height=11, width=3)
survey_graph <- ggplot(melted_survey, aes(Hypothesis, fill=Answer)) +		## Plotted by trait
  geom_bar(position="fill") +
  scale_fill_manual(values=varit) +
  labs(x="", y="Fraction", title="Support for different hypotheses") +
  theme_grey(base_size=9) +
  theme(axis.text.x=element_blank(), plot.title = element_text(hjust = 1.3),
        legend.position = "bottom", legend.direction = "vertical") +
  coord_flip() +
  scale_y_reverse() +
  facet_wrap(~Trait, ncol=1, scales="free")
survey_graph
dev.off()


for(i in unique(melted_survey$Hypothesis)) {							# All hypotheses separately plotted
  temp <- ggplot(melted_survey[melted_survey$Hypothesis == i , ], aes(Answer)) + 
    geom_histogram(fill="green") +
    labs(x = "", y = "count", title = paste(i, " in survey")) +
    theme_grey(base_size=20) +
    coord_cartesian(xlim = NULL, ylim = c(0, 650))
  
  print(temp)
}

#### HISTOGRAMS BY EXPERTISE

deadthing_experts <- ggplot(melted_survey[melted_survey$Expertise=="dead",], aes(Hypothesis, fill=Answer)) +	## Dead thing experts plotted by trait
  geom_histogram() +
  labs(x="", y="Count") +
  labs(title="Dead things experts' views") +
  theme_grey(base_size=24) +
  theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) + 
  facet_wrap(~Result, ncol=4, scales="free_x")

livething_experts <- ggplot(melted_survey[melted_survey$Expertise=="live",], aes(Hypothesis, fill=Answer)) +	## Livething experts plotted by trait
  geom_histogram() +
  labs(x="", y="Count") +
  labs(title="Live things experts' views") +
  theme_grey(base_size=24) +
  theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) + 
  facet_wrap(~Result, ncol=4, scales="free_x")

for(i in unique(melted_survey$Expertise.specific)) {		
  expert <- ggplot(melted_survey[melted_survey$Expertise.specific== i,],		## All specific expertise plotted
                   aes(Hypothesis, fill=Answer)) +
    geom_histogram() +
    labs(x="", y="Count") +
    labs(title=paste(i,"s' views")) +
    theme_grey(base_size=24) +
    theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) +
    facet_wrap(~Trait, ncol=4, scales="free_x")
  
  print(expert)
  
}



for(i in unique(melted_survey$Trait)) {	
  trait <- ggplot(melted_survey[melted_survey$Trait == i,],		## All traits plotted, filled by expertise
                  aes(Answer, fill=Expertise.specific)) +
    geom_histogram() +
    labs(x="", y="Count") +
    labs(title=paste(i,"in survey")) +
    theme_grey(base_size=24) +
    theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) +
    facet_wrap(~Hypothesis)
  
  print(trait)
  
}



for(i in unique(melted_survey$Trait)) {	
  gen_expert <- ggplot(melted_survey[melted_survey$Trait == i,],		## Traits by general expertise
                       aes(Answer, fill=Expertise)) +
    geom_histogram() +
    labs(x="", y="Count") +
    labs(title=paste(i,"in survey")) +
    theme_grey(base_size=24) +
    theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.4)) +
    facet_wrap(~Hypothesis)
  
  print(gen_expert)
  
}


for(i in unique(melted_survey$Trait)) {				## Traits wrapped prob density by specific expertise
  temp <- ggplot(melted_survey[melted_survey$Trait == i , ],
                 aes(x = as.numeric(Answer), colour = Expertise.specific)) + 
    geom_density(adjust = 3) +
    labs(x = "1 = very likely", y = "prob density", title = paste(i, " in survey")) +
    theme_grey(base_size=24) +
    facet_wrap(~ Hypothesis)
  
  print(temp)
}

for(i in unique(melted_survey$Hypothesis)) {			## All hypotheses prob density by specific expertise
  temb <- ggplot(melted_survey[melted_survey$Hypothesis == i , ],
                 aes(x = as.numeric(Answer), colour = Expertise.specific)) + 
    geom_density(adjust = 3) +
    labs(x = "1 = very likely", y = "prob density", title = paste(i, " by expertise")) +
    theme_grey(base_size=24)
  
  print(temb)
}


for(i in unique(melted_survey$Trait)) {				## Traits wrapped prob density by familiarity
  tema <- ggplot(melted_survey[melted_survey$Trait == i , ],
                 aes(x = as.numeric(Answer), colour = Familiarity)) + 
    geom_density(adjust = 3) +
    labs(x = "1 = very likely", y = "prob density", title = paste(i, " in survey")) +
    theme_grey(base_size=24) +
    facet_wrap(~ Hypothesis)
  
  print(tema)
}

}#if(FALSE)

biped.cols <- c("darkorange", "bisque2", "coral1", "brown", "darkgoldenrod1", "firebrick1", "gold4", "turquoise", "indianred1", "tan3")
brain.cols <- c("darkorange", "bisque2", "coral1", "brown", "darkgoldenrod1", "firebrick1", "turquoise", "gold4", "indianred1", "tan3")
hair.cols <- c("coral1", "brown", "darkgoldenrod1", "firebrick1", "turquoise", "gold4", "indianred1", "tan3")

#pdf(file=filename, height=4*mfrow[1], width=4*mfrow[2])
par(mfrow=c(1,4), mar=c(3,3,2,1), mgp=c(1.6,0.6,0))

ggplot(melted_survey[melted_survey$Trait == "Other" , ],			# Other water traits
                 aes(x = as.numeric(Answer), colour = Hypothesis)) +
    geom_freqpoly(binwidth=1, size = 1.5) +
	coord_cartesian(xlim= c(1.5,5.7)) +
	scale_x_discrete(labels = "") +
    labs(x = "very likely on left", y = "count", title = "Other") +
    theme_grey(base_size=24) +
    facet_wrap(~ Expertise, scales="free_y")
#facet_wrap(, ncol=1)

ggplot(melted_survey[melted_survey$Trait == "Bipedalism" , ],			#Bipedalism hypotheses
                 aes(x = as.numeric(Answer), colour = Hypothesis)) +
	scale_colour_manual(values=biped.cols) + 
    geom_freqpoly(binwidth=1, size = 1.5) +
	coord_cartesian(xlim= c(1.5,5.7)) +
	scale_x_discrete(labels = "") +
    labs(x = "very likely on left", y = "count", title = "Bipedalism") +
    theme_grey(base_size=24) +
    facet_wrap(~ Expertise, scales="free_y")

ggplot(melted_survey[melted_survey$Trait == "Brain" , ], 			# Brain hypotheses
                 aes(x = as.numeric(Answer), colour = Hypothesis)) + 
	scale_colour_manual(values=brain.cols) + 
    geom_freqpoly(binwidth=1, size = 1.5) +
	coord_cartesian(xlim= c(1.5,5.7)) +
	scale_x_discrete(labels = "") +
    labs(x = "very likely on left", y = "count", title = "Brain") +
    theme_grey(base_size=24) +
    facet_wrap(~ Expertise, scales="free_y")

ggplot(melted_survey[melted_survey$Trait == "Hairlessness" , ],			# Hairlessness hypotheses
                 aes(x = as.numeric(Answer), colour = Hypothesis)) + 
	scale_colour_manual(values=hair.cols) + 
    geom_freqpoly(binwidth=1, size = 1.5) +
	coord_cartesian(xlim= c(1.5,5.7)) +
	scale_x_discrete(labels = "") +
    labs(x = "very likely on left", y = "count", title = "Hairlessness") +
    theme_grey(base_size=24) +
    facet_wrap(~ Expertise, scales="free_y")

for(i in unique(melted_survey$Trait)) {				## Traits wrapped prob density by familiarity
  tema <- ggplot(melted_survey[melted_survey$Trait == i , ],
                 aes(x = as.numeric(Answer), colour = Hypothesis)) + 
    geom_freqpoly(binwidth=1, size = 1.5) +
	coord_cartesian(xlim= c(1.5,5.7)) +
	scale_x_discrete(labels = "") +
    labs(x = "very likely on left", y = "count", title = paste(i)) +
    theme_grey(base_size=24) +
    facet_wrap(~ Expertise, scales="free_y")
  
  print(tema)
}

Principal coordinates and component analyses, linear regression

Model run [1] with linear regression
Important functions: princomp, biplot.princomp

+ Show code - Hide code

# This is code Op_en7839/ on page [[How scientists perceive the evolutionary origin of human traits: results of a survey study]]

library(OpasnetUtils)
library(ggplot2)
library(reshape2)
library(vegan)

objects.latest("Op_en7839", code_name="publishdata") # This line downloads data.frames with the data

hypocols1 <- 29:79

# Colname changes should happen already in the the beginning in the t2b.

colnames(survey1) <- gsub("^[0-9]*\\.", "", colnames(survey1))

survey <- survey1[survey1$Allhypoanswers , ]

answ <- as.matrix(survey[hypocols1])

if(FALSE){

########### Linear regression

#survey1$Occupation <- factor(survey1$Occupation, ordered = FALSE)

# Multiple Linear Regression Example
fit <- lm(aah.average ~ Gender + Age + Education + Pub.all.peer + Pub.hum.peer + Pub.hum.popular, data = survey1)
fit <- lm(aah.average ~ Gender + Teaching + Familiarity + Expertise.specific , data = survey1)
fit <- lm(aah.average ~ Gender + Expertise , data = survey1)
#fit <- lm(aah.average ~ Gender + Occupation , data = survey1)
fit <- lm(aah.average ~ Experience + Education , data = survey1)
# The last analysis contains all interesting variables.
fit <- lm(aah.average ~ Gender + Pub.all.peer + Expertise.specific + Familiarity + Teaching , data = survey1)

survey1$Highpublications <- survey1$Pub.all.peer >= "11-40"
fit <- lm(aah.average ~ Gender + Highpublications + Expertise.specific + Familiarity + Teaching , data = survey1)
oprint(summary(fit)) # show results
# Important variables that increase acceptability: Gender (females), Expertise.specific (biologist, 
# geologist, human biologist, other), familiarity (less), teaching (more)

tapply(survey1$Pub.all.peer, survey1$Pub.all.peer, FUN = length)
# Other useful functions
oprint(coefficients(fit)) # model coefficients
oprint(confint(fit, level=0.95)) # CIs for model parameters
#fitted(fit) # predicted values
#residuals(fit) # residuals
oprint(anova(fit)) # anova table
#vcov(fit) # covariance matrix for model parameters
#influence(fit) # regression diagnostics 

survey1$all.average <- 5 - rowSums(survey1[c(
  "Sociality", 
  "Language",
  "Articulation",
  "Energy.supply",
  "Tool.use",
  "Social"
)])/6

# Multiple Linear Regression Example
#fit <- lm(all.average ~ Gender + Age + Education + Pub.all.peer + Pub.hum.peer + Pub.hum.popular, data = survey1)
#fit <- lm(all.average ~ Gender + Teaching + Familiarity + Expertise.specific , data = survey1)
#fit <- lm(all.average ~ Gender + Expertise , data = survey1)
#fit <- lm(all.average ~ Gender + Occupation , data = survey1)
#fit <- lm(all.average ~ Experience + Education , data = survey1)
#fit <- lm(all.average ~ Gender + Pub.all.peer + Expertise.specific + Familiarity + Teaching , data = survey1)
# The last analysis contains all interesting variables.
fit <- lm(all.average ~ Gender + Highpublications + Expertise.specific + Familiarity + Teaching , data = survey1)
oprint(summary(fit)) # show results

survey1$Allsupport <- 5 - rowSums(survey1[29:79]) / 51
fit <- lm(Allsupport ~ Expertise.specific , data = survey1)

survey1$nonaah.average <- 5 - rowSums(survey1[c(29:30, 32:52, 54:57, 59:62, 64, 66:70)])/37 # Aah: 31 53 58 63 65 71-]) / 51
fit <- lm(nonaah.average ~ Expertise.specific , data = survey1)
oprint(summary(fit))
nonaah.averagefit <- lm(aah.average ~ Expertise.specific , data = survey1)
oprint(summary(nonaah.averagefit))

test <- as.data.frame(lapply(survey1[5:19], FUN = function(x) !is.na(x)))
test2 <- rowSums(test) 
survey1[5:19] <- as.data.frame(lapply(test, FUN = function(x) ifelse(test2 == 0, NA, x)))

fit <- lm(aah.average ~ Anthropol + Bio.physiol + Bio.ecol + Bio.evol + Bio.genet + Bio.other +
            Geology + Hum.cardiov + Hum.muscle + Hum.neurol + Hum.nutr + Hum.other +
            Paleoanthr + Paleontol + Other.open , data = survey1)
oprint(summary(fit))


### PRINCIPAL COORDINATES ANALYSIS

# run PCoA with people as the units of observation
pcoa_people <- cmdscale(dist(answ))
pcoa_people <- as.data.frame(pcoa_people)

pcoa_people$familiarity <- survey$Familiarity
pcoa_people$Expertise <- survey$Expertise.specific
pcoa_people$aah.average <- survey$aah.average
pcoa_people$publication <- survey$Pub.all.peer
pcoa_people$he.publication <- survey$Pub.hum.peer


for(i in unique(pcoa_people$Expertise)) {
	kuva <- ggplot(pcoa_people, aes(x=V1, y=V2, colour=familiarity)) + 		 # every dot is a person
	geom_point(aes(size=ifelse(pcoa_people$Expertise == i, 1, 0))) +
	scale_size(range=c(1,5)) +
	theme_grey(base_size=24)+
	labs(x="", y="", title=paste(i))

	print(kuva)
}

ggplot(pcoa_people, aes(x=V1, y=V2, colour=Expertise)) + 		 # every dot is a person
	geom_point(aes(size=familiarity)) +
	theme_grey(base_size=24)+
	labs(x="", y="")


ggplot(pcoa_people, aes(x=V1, y=V2)) + 		 # every dot is a person, aah opinion sze
	geom_point(aes(size=aah.average)) +
	scale_size(range=c(1,5)) +
	theme_grey(base_size=24)+
	labs(x="", y="")


ggplot(pcoa_people, aes(x=V1, y=V2, color=publication)) + 		 # every dot is a person, publications size and color
	geom_point(aes(size=publication)) +
	theme_grey(base_size=24)+
	labs(x="", y="")

ggplot(pcoa_people, aes(x=V1, y=V2, color=publication)) + 		 # every dot is a person, publications color
	geom_point(size=3) +
	theme_grey(base_size=24)+
	labs(x="", y="")

ggplot(pcoa_people, aes(x=V1, y=V2)) + 		 # every dot is a person, human evolution publications size
	geom_point(aes(size=he.publication)) +
	theme_grey(base_size=24)+
	labs(x="", y="")

} #if(FALSE)

# run PCoA with hypotheses as the units of observation
pcoa_hypot <- cmdscale(dist(t(answ)))
pcoa_hypot <- as.data.frame(pcoa_hypot)
pcoa_hypot$trait <- c(rep("bipedalism", 10), rep("brain", 10), rep("hairlessness", 8), rep("fat", 4), rep("larynx", 3), rep("speech", 7), rep("other", 9))
pcoa_hypot$hypothesis <- rownames(pcoa_hypot)

pcoa_hypot$scale <-  5 - colMeans(answ)

ggplot(pcoa_hypot, aes(x=V1, y=V2, colour = trait, label = hypothesis)) +  # every dot is a hypothesis
	geom_point(aes(size=scale)) +
	scale_size(range=c(3,15)) +
	theme_grey(base_size=24) +
	labs(x="", y="") +
	geom_text(check_overlap = TRUE, colour="black")

#rda_some <- rda(t(answ))
#plot(rda_some)

## PCoA with capscale

aah.pop.sizes <- survey$aah.average4/2
pop12.sizes <- survey$nonaah.average12/2
pop6.sizes <- survey$nonaah.average6/2
social.sizes <- survey$Social.average/2
enviro.sizes <- survey$Enviro.average/2

answ2 <- 5-answ

pcoa_caps <- capscale(t(answ2) ~ 1, distance="euclidean") ##PCoA done

sizes <- survey$aah.average/2
traits <- as.factor(pcoa_hypot$trait)
colstr <- c("palevioletred1","royalblue1","seagreen1","violet","khaki2","skyblue", "orange")
trait.cols <- colstr[traits]

plot(pcoa_caps, display = c("sp", "wa"), type="n") ## PCoA biplot, full scale
points(pcoa_caps, display= c("sp"), col="gray40") # adding the people points
points(pcoa_caps, display= c("wa"), pch=19, cex=(pcoa_hypot$scale*1.7), col=trait.cols)
text(pcoa_caps, display=c("wa"), srt=25, cex=0.7)
legend("topleft", levels(traits), fill=colstr, bty="n")

familiarity <- as.factor(survey$Familiarity)
publications <- as.factor(survey$Pub.all.peer)
expertise <- as.factor(survey$Expertise)
gender <- as.factor(survey$Gender)
colspub <- c("seagreen1","khaki2","palevioletred1","royalblue1","black")
colsfml <- c("khaki2","palevioletred1","royalblue1")
colsexp <- c("seagreen1","khaki2","palevioletred1","royalblue1")
colsgen <- c("palevioletred1","royalblue1")
pub.colors <- colspub[publications]
famil.colors <- colsfml[familiarity]
expert.colors <- colsexp[expertise]
gender.colors <- colsgen[gender]

plot(pcoa_caps, display = c("wa"), xlim=c(-0.7,0.5), ylim = c(-0.4,0.6)) ## Biplot people, aah average sized, publications color
points(pcoa_caps, display = c("sp"), pch=19, cex=sizes, col=pub.colors)
legend(0.25,0.4, levels(publications), fill=colspub, bty="n")

plot(pcoa_caps, display = c("wa"), xlim=c(-0.7,0.5), ylim = c(-0.4,0.6)) ## Biplot people, aah average sized, falimiarity color
points(pcoa_caps, display = c("sp"), pch=19, cex=sizes, col=famil.colors)
legend(0.25,0.4, levels(familiarity), fill=colsfml, bty="n")

plot(pcoa_caps, display = c("wa"), xlim=c(-0.7,0.5), ylim = c(-0.4,0.6)) ## Biplot people, aah average sized, expertise color
points(pcoa_caps, display = c("sp"), pch=19, cex=sizes, col=expert.colors)
legend(0.25,0.4, levels(expertise), fill=colsexp, bty="n")

plot(pcoa_caps, display = c("wa"), xlim=c(-0.7,0.5), ylim = c(-0.4,0.6)) ## Biplot people, aah average sized, expertise color
points(pcoa_caps, display = c("sp"), pch=19, cex=sizes, col=gender.colors)
legend(0.25,0.4, levels(gender), fill=colsgen, bty="n")

plot(pcoa_caps, display = c("wa"),					## Biplot people, aah average sized, circles
	xlim=c(-0.7,0.5), ylim = c(-0.4,0.6), main="AAH-average as size")
points(pcoa_caps, display = c("sp"), pch=21, cex=sizes)

plot(pcoa_caps, display = c("wa"),					## Biplot people, aah popular sized, circles
	xlim=c(-0.7,0.5), ylim = c(-0.4,0.6), main="4 most popular AAH average as size")
points(pcoa_caps, display = c("sp"), pch=21, cex=aah.pop.sizes)

plot(pcoa_caps, display = c("wa"),					## Biplot people, non-aah 12 popular sized, circles
	xlim=c(-0.7,0.5), ylim = c(-0.4,0.6), main="12 most popular non-AAH average as size")
points(pcoa_caps, display = c("sp"), pch=21, cex=pop12.sizes)

plot(pcoa_caps, display = c("wa"),					## Biplot people, non-aah 6 popular sized, circles
	xlim=c(-0.7,0.5), ylim = c(-0.4,0.6), main="6 most popular non-AAH average as size")
points(pcoa_caps, display = c("sp"), pch=21, cex=pop6.sizes)

plot(pcoa_caps, display = c("wa"),					## Biplot people, non-aah social-related sized, circles
	xlim=c(-0.7,0.5), ylim = c(-0.4,0.6), main="Non-AAH social-related average as size")
points(pcoa_caps, display = c("sp"), pch=21, cex=social.sizes)

plot(pcoa_caps, display = c("wa"),					## Biplot people, no-aah environment-related sized
	xlim=c(-0.7,0.5), ylim = c(-0.4,0.6), main="Non-AAH environment-related average as size")
points(pcoa_caps, display = c("sp"), pch=21, cex=enviro.sizes)

if(FALSE){

###### Principal component analysis

prin <- princomp(answ, cor = FALSE)

plot(prin)

biplot(prin, display = c("sp", "wa"))

} #if(FALSE)

AAH criticism and anova

Model run 8.9.2017 with all analyses [2]

The code contins the following parts in this order:

Anova and Tukey's test
Support and rejection of AAH in literature
Number of credible hypotheses per trait (individual respondent's view)
Explanations for not finding any hypothesis likely
Means, SD, and support distributions for all hypotheses ~ Expertise
Support distribution for all AAH critiques
Credibility of the two most likely hypothesis on the respondent level
Support distribution for the most and second moset likely hypotheses
Mean and SD of average indices by different determinants
Distribution of credibility scores by Expertise
AAH criticism
Hypothesis conflicts (based on the t2b conflict tables)
Number of credible hypotheses
Attitudes observed in the literature
PCoA of criticism
Criticism in AAH familiarity subgroups
Mean and SD of individual hypotheses by different determinants

+ Show code - Hide code

# This is code Op_en7839/ on page [[How scientists perceive the evolutionary origin of human traits: results of a survey study]]

library(OpasnetUtils)
library(ggplot2)
library(reshape2)
library(vegan)

# Colour palettes for graphs
pal5 <- c("seagreen1","khaki2","darkorange1","darkblue","black")
pal3 <- c("seagreen1","darkorange1","darkblue")
pal4 <- c("#8DD3C7", "#FFFFB3", "#BEBADA", "#FB8072") # From ColorBrewer Set3 (as in Fig 5)
# making colour palette for values 1.0, 1.1, ... , 5.0
pal41 <- colorRampPalette(c("indianred1", "royalblue3"))(41) # 41-step colour slider
varit4 <- c("mediumseagreen", "greenyellow", "grey70", "darkgoldenrod1", "darkorange")

objects.latest("Op_en7839", code_name="publishdata") # This line downloads data.frames with the data

colnames(survey1)

############################################ How many have heard before about AAH

table(survey1$Hear.before)

############################################# Larynx and Articulation equally likely

sum(survey1$Articulation == survey1$Larynx.S & as.numeric(survey1$Articulation) <= 2, na.rm = TRUE)/nrow(survey1)

################################# ANOVA AND TUKEY'S TEST

#### test if means differ among gender, age or expertise within aah_average and nonaah_average and make a table

# make a list to store statistical results
varnames <- c("aah_average", "nonaah_average", "nonaah.average12", "Social.average", "Enviro.average")
explnames <- c("Gender", "Age", "Expertise", "Expertise2", "Pub.all.peer", "Pub.hum.peer", "Pub.all.popular",
               "Pub.hum.popular", "Familiarity", "Teaching")
M <- matrix(data=NA, nrow=length(explnames), ncol=length(varnames))
colnames(M) <- varnames
rownames(M) <- explnames
res.all <- apply(M, 1, as.list)

# run the analyses of variance
for(i in explnames) {
  for(j in varnames) {
    subset <- survey1[,c(i,j)]
    colnames(subset) <- c("Explanation", "Variable")
    #    subset$Variable <- as.numeric(subset$Variable)
    res.aov <- aov(Variable ~ Explanation, data=subset)
    res.tukey <- TukeyHSD(res.aov)
    res.all[[i]][[j]][1:2] <- list(res.aov, res.tukey)
  }
} 

oprint(res.all)

####################### Means and SDs of average indices by several determinants

cat("Means and SDs of average indices by several determinants.\n")

meantable <- function(col) {
  out <-  aggregate(
    survey1[varnames], 
    by = survey1[col], 
    FUN = function(x) {
        paste(
          sprintf("%.2f", mean(as.numeric(x), na.rm = TRUE)), " (",
          sprintf("%.2f", sd(as.numeric(x), na.rm = TRUE)), ")", sep=""
          )
    }
  )

  colnames(out)[1] <- "Subgroup"
  out <- rbind(
    data.frame( # significant p values summarised
      Subgroup = col, 
        sapply(
        varnames,
        FUN = function(x) {
          pees <- res.all[[col]][[x]][2][[1]][[1]][,4] # p values from ANOVA/Tukey
          pees <- pees[pees < 0.05]
          pees <- paste(
            paste(
              names(pees),
              cut(pees, breaks = c(-1,0.001,0.01, 0.05), labels = c("***", "**", "*"))
            ), collapse = " "
          )
        },
        simplify = FALSE
      )
    ),
    out # means and sd's
  )
  return(out)
}
means <- rbind(
  meantable("Expertise"),
  meantable("Familiarity"),
  meantable("Gender"),
  meantable("Age"),
  meantable("Pub.hum.peer"),
  meantable("Teaching")
)

oprint(means)

############################################# ANOVA TUKEY WITH INTERACTIONS

res.aov <- aov(aah_average ~ Teaching*Expertise, data=survey1)
res.tukey <- TukeyHSD(res.aov)
oprint(list(res.aov, res.tukey))

res.aov <- aov(aah_average ~ Familiarity*Expertise, data=survey1)
res.tukey <- TukeyHSD(res.aov)
oprint(list(res.aov, res.tukey))

############################################# TUKEY INTERACTIONS IN NICE TABLES

res.aov <- aov(aah_average ~ Familiarity*Expertise, data=survey1)
res.tukey <- TukeyHSD(res.aov)
out <- res.tukey[[3]][order(dimnames(res.tukey[[3]])[[1]]),]

test <- lapply( # Split the comparison items
  strsplit(dimnames(out)[[1]], split="-"), 
  FUN = function(x) strsplit(x, split = ":")
)
oprint(out[sapply(test, FUN = function(x) x[[1]][1] == x[[2]][1]),]) # Matching Familiarity
oprint(out[sapply(test, FUN = function(x) x[[1]][2] == x[[2]][2]),]) # Matching Expertise


############################# Attitudes observed in the literature

ggplot(survey1, aes(x=Expertise, fill=Rejecting))+geom_bar(position="stack")+coord_flip()
ggplot(survey1, aes(x=Expertise, fill=Neutral))+geom_bar(position="stack")+coord_flip()
ggplot(survey1, aes(x=Expertise, fill=Supportive))+geom_bar(position="stack")+coord_flip()

ggplot(survey1, aes(x=Familiarity, fill=Rejecting))+geom_bar(position="stack")+coord_flip()
ggplot(survey1, aes(x=Familiarity, fill=Neutral))+geom_bar(position="stack")+coord_flip()
ggplot(survey1, aes(x=Familiarity, fill=Supportive))+geom_bar(position="stack")+coord_flip()

oprint(table(survey1[c("Familiarity", "Rejecting")], useNA = "ifany"))
oprint(table(survey1[c("Familiarity", "Neutral")], useNA = "ifany"))
oprint(table(survey1[c("Familiarity", "Supportive")], useNA = "ifany"))

ggplot(survey1, aes(x=aah_average, colour=Rejecting))+geom_density()
ggplot(survey1, aes(x=aah_average, colour=Neutral))+geom_density()
ggplot(survey1, aes(x=aah_average, colour=Supportive))+geom_density()

########################### Number of credible hypotheses

hypos <- allhypos[-(71:79-28)] #"Other" trait removed because not meaningful here.
traits <- as.character(survey1_questions$Trait[
  match(
    hypos, 
    survey1_questions$Name
  )
])

numHypotheses <- function(credibility) {
  # support is an array for tests whether a trait*respondent pair
  # fulfils the given credibility criteria (1:very likely, 5: very unlikely).
  support <- sapply(survey1[hypos], FUN = as.numeric) <= 6-credibility
  colnames(support) <- gsub("\\.", " ", colnames(support))
  
  # test is an array with number of credible hypothesis for each 
  # respondent * trait pair.
  test <- sapply(
    unique(traits), FUN = function(x) {
      rowSums(support[,traits %in% x], na.rm = TRUE)
    }
  )
  # Combine with the determinant columns in survey1 and melt.
  test <- cbind(survey1[c(1:4,23:28,123)], test)
  test <- melt(
    test, 
    measure.vars = unique(traits), 
    variable.name = "Trait", 
    value.name = "Credible.hypotheses"
  )
  return(test)
}

numhypos <- rbind(
  data.frame(numHypotheses(5), Likelihood = "Very likely"),
  data.frame(numHypotheses(4), Likelihood = "Likely")
)
numhypos$Credible.hypotheses <- factor(
  numhypos$Credible.hypotheses, 
  levels = as.character(0:10),
  ordered = TRUE
)

numhypos <- merge(
  numhypos, 
  aggregate(survey1_questions$Name, by = survey1_questions["Trait"], length)
)
numhypos$Trait2 <- paste(numhypos$Trait," (",numhypos$x,")", sep="")

oprint(
  reshape(
    aggregate(
      numhypos$Credible.hypotheses, 
      by = numhypos[c("Credible.hypotheses", "Trait2", "Likelihood")], 
      function(x) length(x)/nrow(survey1)
    ),
    timevar = "Likelihood",
    idvar = c("Credible.hypotheses", "Trait2"),
    v.names = "x",
    direction = "wide"
  )
)

ggplot(numhypos, aes(x=Trait2, fill=Credible.hypotheses))+
  geom_bar(position="fill")+coord_flip()+ facet_wrap(~ Likelihood)+
  scale_fill_brewer(palette = "Set3", name="Number of\nhypotheses")+
  labs(title = "Counting likely or very likely hypotheses",
       y = "Percentage of respondents")+
  theme(axis.text.y = element_text(size=10), 
        legend.text= element_text(size=10),
        axis.text.x = element_blank())+
  scale_y_reverse()

################### Explanations for not finding any hypothesis likely

oprint(
    data.frame(`Explanations for not believing any hypotheses` = unique(
    numhypos[
      numhypos$Likelihood == "Likely" & numhypos$Credible.hypotheses == 0 , 
      "Other.comments"
    ]
  ))
)

###################### Means, SD, and support distributions for all hypotheses

cat("Means and SD of support for each hypothesis. 1 = very unlikely, 5 = very likely.\n")

out <- aggregate(
  melted_survey["Answer"], 
  by = melted_survey[c("Hypothesis", "Trait")], 
  FUN = function(x) c(
    Mean = mean(6 - as.numeric(x), na.rm = TRUE),
    SD = sd(as.numeric(x), na.rm = TRUE)
  )
)

out <- cbind(out[1:2], out[[3]])
oprint(out)

cat("Support distribution for all hypotheses.\n")

frek <- aggregate(
  melted_survey$Answer, 
  by = melted_survey[c("Answer", "Hypothesis", "Trait")], 
  FUN = length
)
frek <- aggregate(
    frek$x, 
    by = frek[c("Hypothesis", "Trait")], 
    FUN = function(x) x/sum(x)
)
colnames(frek[[3]]) <- levels(melted_survey$Answer)
frek <- cbind(frek[1:2], frek[[3]])
oprint(frek)

################################################ Support distribution for all AAH critiques

cat("Support distribution for all AAH critiques.\n")

oprint(table(critique[c("Critique", "Answer")], useNA = "ifany"))

oprint(
  tapply(
    critique$Critique, 
    INDEX = critique[c("Critique", "Answer")],
    FUN = function(x) length(x)/nrow(survey1)
  )
)

########################## Credibility of the two most likely hypothesis on the respondent level

frek <- melted_survey
frek$Trait2 <- ifelse(
  frek$Trait == "Other", 
  as.character(frek$Hypothesis), 
  as.character(frek$Trait)
)

frek.1 <- aggregate(
  6-as.numeric(frek$Answer), 
  by = frek[c("ID", "Trait2", "Expertise", "Gender", "Age", "Familiarity")], 
  FUN = max
)

frek.2 <- aggregate(
  6-as.numeric(frek$Answer), 
  by = frek[c("ID", "Trait2", "Expertise", "Gender", "Age", "Familiarity")], 
  FUN = function(x) -sort(-x)[2]
)
freks <- rbind(
  data.frame(Hypothesis_rank = "Best", frek.1),
  data.frame(Hypothesis_rank = "Second best", frek.2)
)
freks <- freks[!is.na(freks$x),]

ggplot(freks, aes(x=as.character(x), fill = Hypothesis_rank))+geom_bar(position = "dodge")+
  facet_wrap(~ Trait2)+
  labs(
    x="Credibility score",
    title="Credibility of the two best hypothesis on respondent level")+
  theme_gray(base_size=24)

################################ Support distribution for the most and second moset likely hypotheses

cat("Support distribution of the most and second most likely hypotheses for all traits.\n")

frek2 <- melt(
  tapply(
    freks$x, 
    INDEX = freks[c("x", "Trait2", "Hypothesis_rank")], 
    FUN = length
  )
)
frek2$value[is.na(frek2$value)] <- 0
frek2 <- aggregate(
  frek2$value, 
  by = frek2[c("Trait2", "Hypothesis_rank")], 
  FUN = function(x) x/sum(x)
)
frek2 <- cbind(frek2[1:2], frek2[[3]])
colnames(frek2)[3:7] <- rev(levels(melted_survey$Answer))
frek2 <- frek2[!is.nan(frek2[[3]]),]

oprint(frek2)

freks3 <- frek[c("Hypothesis", "ID", "Trait2", "Answer")]
freks3$x <- 6-as.numeric(freks3$Answer)
freks3$IDm <- 1:nrow(freks3)

cat("Frequencies of pairs of two most credible hypotheses (on respondent level) for each trait.\n")

freks5 <- aggregate(
  freks3$IDm,
  by = freks3[c("ID", "Trait2")],
  FUN = function(x) {
    temp <- freks3[x , c("Hypothesis", "x")]
    paste(sort(temp$Hypothesis[order(-temp$x)[1:2]]), collapse = " - ")
  }
)
freks5 <- aggregate(freks5$x, by = freks5[c("Trait2", "x")], FUN = length)
freks5 <- freks5[order(freks5[[1]], -freks5[[3]]) , ]

oprint(freks5)

####################################### Distribution of credibility scores by Expertise

cat("Support distribution for all hypotheses by expertise.\n")

frek <- melt(
  tapply(
    melted_survey$Answer, 
    INDEX = melted_survey[c("Answer", "Hypothesis", "Trait", "Expertise")], 
    FUN = length
  )
)
frek$value[is.na(frek$value)] <- 0
frek <- aggregate(
  frek$value, 
  by = frek[c("Hypothesis", "Trait", "Expertise")], 
  FUN = function(x) x/sum(x)
)
frek <- cbind(frek[1:3], frek[[4]])
colnames(frek)[4:8] <- levels(melted_survey$Answer)
frek <- frek[!is.nan(frek[[4]]),]

oprint(frek)

######################## AAH criticism 

# AAH criticism by expertise
ggplot(critique, aes(x = Critique, fill = Answer))+
  geom_bar(position="fill")+coord_flip()+
  facet_grid(. ~ Expertise)+theme_gray(base_size = 24)+theme(legend.position = "bottom")+
  scale_fill_manual(values=varit4)


#AAH criticism in the whole group

ggplot(critique, aes(x = Critique, fill = Answer))+
  geom_bar(position="fill")+coord_flip()+
  labs(x = "", y = "Fraction")+
  theme_gray(base_size = 24)+
  scale_fill_manual(values=varit4)

##############AAH criticism in the whole group with Heard-before blocks

ggplot(critique, aes(x = Critique, fill=Answer, color=Hear.before)) +	
  geom_bar(position="fill") +
  scale_fill_manual(values=varit4) +
  scale_color_manual(values = c("navyblue", "#FFFFFF00")) +
  labs(x="", y="Fraction", title="Support for AAH criticism") +
  theme_grey(base_size=9) +
  coord_flip() +
  scale_y_reverse()

####################### Hypothesis conflicts

# function conflict.test produces, for each trait and respondent,
# an indicator about whether there are two conflicting hypotheses at a given
# credibility level or higher. 
# Whether two hypotheses are conflicting or whether they
# can coincide, is based on authors' evaluation and this data is stored
# in t2b tables with the name of each trait. If a conflict is not clear, 
# it is marked with "?" in the tables. If sensitivity to conflict is low,
# "?" equals "n", and if it is high, it equals "y".
# These are in different columns in data.frame conflict.
# The function takes the name of column as argument "Sensitivity"

traits <- as.character(survey1_questions$Trait[
   match(
      conflict$Hypothesis1, 
      gsub("\\.", " ", survey1_questions$Name)
   )
])

conflict.test <- function(credibility, sensitivity) {
  # support is an array for tests whether a trait*respondent pair
  # fulfils the given credibility criteria.
  support <- sapply(survey1[allhypos], FUN = as.numeric) >= credibility
  colnames(support) <- gsub("\\.", " ", colnames(support))

  # test is an array with conflict test for each 
  # respondent * pair of traits in conflict.
  test <- sapply(
    1:nrow(conflict), FUN = function(x) {
      support[,conflict$Hypothesis1[x]] &
      support[,conflict$Hypothesis2[x]] &
      conflict[[sensitivity]][x]
    }
  )
  # aggregate over pairs to the traits (assuming that both hypotheses explain
  # the same trait)
  test <- aggregate(
    data.frame(t(test)), 
    by = list(traits),
    FUN = any
  )
  tr <- test[[1]]
  test <- t(test[-1])
  colnames(test) <- tr
  # Combine with the determinant columns in survey1 and melt.
  test <- cbind(survey1[1:28], test)
  test <- melt(
    test, 
    measure.vars = tr, 
    variable.name = "Trait", 
    value.name = "Conflict"
  )
  return(test)
}

test <- conflict.test(credibility = 4, sensitivity = "Conflict.hi")

ggplot(test, aes(x = Trait, fill = Conflict))+geom_bar(position="fill")+
  coord_flip()+
  labs(title="Conflict between two credible hypotheses", y = "Fraction")

test <- conflict.test(credibility = 4, sensitivity = "Conflict.lo")

ggplot(test, aes(x = Trait, fill = Conflict))+geom_bar(position="fill")+
  coord_flip()+
  labs(title="Conflict between two credible hypotheses", y = "Fraction")

test <- conflict.test(credibility = 5, sensitivity = "Conflict.hi")

ggplot(test, aes(x = Trait, fill = Conflict))+geom_bar(position="fill")+
  coord_flip()+
  labs(title="Conflict between two very credible hypotheses", y = "Fraction")

test <- conflict.test(credibility = 5, sensitivity = "Conflict.lo")

ggplot(test, aes(x = Trait, fill = Conflict))+geom_bar(position="fill")+
  coord_flip()+
  labs(title="Conflict between two credible hypotheses", y = "Fraction")

############################## Principal coordinates analysis on hypotheses

answ <- 6 - sapply(survey1[survey1$Allhypoanswers,allhypos], as.numeric)

pcoa_caps <- capscale(t(answ) ~ 1, distance="euclidean") ##PCoA done

familiarity <- as.factor(survey1$Familiarity)
publications <- as.factor(survey1$Pub.all.peer)
expertise <- as.factor(survey1$Expertise)
aah.average <- as.factor(survey1$aah_average)
colspub <- c("seagreen1","khaki2","darkorange1","darkorchid","black")
colsfml <- c("seagreen1","darkorange1","darkorchid")
colsexp <- c("seagreen1","khaki2","darkorange1","darkorchid")
colaah <- colorRampPalette(c("indianred1", "royalblue3")) # making the gradient for colors
pub.colors <- colspub[publications]
famil.colors <- colsfml[familiarity]
expert.colors <- colsexp[expertise]
col.aah <- colaah(53)[aah.average]
aah_average <- as.numeric(aah.average)

aah.legend <- c(min(aah_average), round(median(aah_average), 2), max(aah_average))
legend.color <- c(col.aah[match(min(aah_average), aah_average)],
                  col.aah[match(median(aah_average), aah_average)],
                  col.aah[match(max(aah_average), aah_average)])

legend.color2 <- c(col.aah[match(5, aah_average)],
                   col.aah[match(4, aah_average)],
                   col.aah[match(3, aah_average)],
                   col.aah[match(2, aah_average)],
                   col.aah[match(1, aah_average)])


plot(pcoa_caps, display = c("wa"), type="n", xlim=c(-0.6,0.5), ylim = c(-0.6,0.5), main="A) Average AAH score") ## Biplot people, publications color
points(pcoa_caps, display = c("sp"), pch=21, col=col.aah)
#legend(0.05,0.5, c("5", "4", "3", "2", "1"), fill=legend.color2, bty="n")
#legend(0.05,0.5, c("Very unlikely", "Unlikely", "No opinion", "Likely", "Very likely"), fill=legend.color2, bty="n")
#legend(0.05,0.5, c("AAH likely", "AAH unlikely"), fill=c("royalblue3", "indianred1"), bty="n")

########### Extra figure 2E about the 12 most popular hypotheses

survey <- survey1[survey1$Allhypoanswers , ]

plot(pcoa_caps, display = c("wa"), type="n", xlim=c(-0.6,0.5), ylim = c(-0.6,0.5), 
     main="E) 12 popular non-AAH hypotheses") ## Biplot people, publications color
points(pcoa_caps, display = c("sp"), pch=21, 
       col=pal41[survey$nonaah.average12*10-9])
legend(0.05,0.5, levels(survey$Wading), fill=pal41[(5:1)*10-9], bty="n")

######################### Principal coordinates analyses of criticism

answ <- 6 - sapply(survey1[89:108], as.numeric)
answ <- answ[!is.na(rowSums(answ)),]

pcoa_caps <- capscale(t(answ) ~ 1, distance="euclidean") ##PCoA done

plot(pcoa_caps, display = c("sp", "wa"), type="n", xlim=c(-6,4.5)) ## PCoA biplot, full scale
points(pcoa_caps, display= c("sp"), col="gray40") # adding the people points
points(pcoa_caps, display= c("wa"), pch=19, cex=colMeans(answ), col="skyblue")
text(pcoa_caps, display=c("wa"), srt=25, cex=1)

plot(pcoa_caps, display = c("wa"), type="n", xlim=c(-0.6,0.5), ylim = c(-0.6,0.5), 
     main="Respondents and criticism") ## Biplot people, publications color
points(pcoa_caps, display = c("sp"), pch=21, col=pal41[rowMeans(answ)*10-9])#, col="navyblue")
legend(-1,0.5, c("5", "4", "3", "2", "1"), fill=pal41[(5:1)*10-9], bty="n")

############################################ FIG 3

ggplot(melted_survey, aes(Hypothesis, fill=Answer)) +		## Plotted by trait
  geom_bar(position="fill") +
  scale_fill_manual(values=varit4) +
  labs(x="", y="Percentage of respondents", title="Support for different hypotheses") +
  theme_grey(base_size=9) +
#  theme(axis.text.x=element_blank(), plot.title = element_text(hjust = 0),
#        legend.position = "bottom", legend.direction = "horizontal") + #  legend.text.align = 0,
#  coord_flip() +
#  scale_y_reverse() + 
  facet_wrap(~Trait, ncol=2, scales="free")

####################### Critique figure in subgroups that did or did not know AAH beforehand

oprint(table(critique[critique$Critique=="Hairy aquatics",1:2]))

ggplot(critique[critique$Hear.before == "Yes",], aes(x = Critique, fill=Answer)) +	
  geom_bar(position="fill") +
  scale_fill_manual(values=varit4) + facet_grid(. ~ Expertise)+
  labs(x="", y="Percentage of respondents", title="Support for AAH criticism (those that had heard of AAH)") + 
  theme_grey(base_size=9) +
  theme(#axis.text.x=element_blank(), #plot.title = element_text(hjust = 1.3), 
        legend.position = "bottom", legend.text.align = 0, legend.direction = "horizontal") +
  coord_flip() #+
#  scale_y_reverse() # This code works on own computer but not in Opasnet. Probably due to a different ggplot2 version

ggplot(critique[critique$Hear.before == "No",], aes(x = Critique, fill=Answer)) +	
  geom_bar(position="fill") +
  scale_fill_manual(values=varit4) + facet_grid(. ~ Expertise)+
  labs(x="", y="Percentage of respondents", title="Support for AAH criticism (those that had not heard of AAH)") +
  theme_grey(base_size=9) +
  theme(#axis.text.x=element_blank(), #plot.title = element_text(hjust = 1.3), 
        legend.position = "bottom", legend.text.align = 0, legend.direction = "horizontal") +
  coord_flip() #+
#  scale_y_reverse() # This code works on own computer but not in Opasnet. Probably due to a different ggplot2 version

#####################################################
# Means and SDs for each hypothesis. Should be merged with the code above.

varnames <- colnames(survey1)[29:79]
explnames <- c("Gender", "Age", "Expertise", "Expertise2", "Pub.all.peer", "Pub.hum.peer", "Pub.all.popular",
               "Pub.hum.popular", "Familiarity", "Teaching")
M <- matrix(data=NA, nrow=length(explnames), ncol=length(varnames))
colnames(M) <- varnames
rownames(M) <- explnames
res.all <- apply(M, 1, as.list)

# run the analyses of variance
for(i in explnames) {
  for(j in varnames) {
    subset <- survey1[,c(i,j)]
    colnames(subset) <- c("Explanation", "Variable")
    subset$Variable <- as.numeric(subset$Variable)
    res.aov <- aov(Variable ~ Explanation, data=subset)
    res.tukey <- TukeyHSD(res.aov)
    res.all[[i]][[j]][1:2] <- list(res.aov, res.tukey)
  }
} 

cat("Means and SDs of individual hypotheses by several determinants.\n")

means <- rbind(
  meantable("Expertise"),
  meantable("Familiarity"),
  meantable("Gender"),
  meantable("Age"),
  meantable("Pub.hum.peer"),
  meantable("Pub.all.peer"),
  meantable("Teaching")
)

oprint(means)

Images

This is the code used to draw the final images. Does not work well in Opasnet at all, but can be copied directly to R.

+ Show code - Hide code

# This is code Op_en7839/ on page [[How scientists perceive the evolutionary origin of human traits: results of a survey study]]

library(ggplot2)
library(OpasnetUtils)
library(reshape2)
library(vegan)
library(gridExtra)

objects.latest("Op_en7839", code_name="publishdata") # This line downloads data.frames with the data

varit <- c("mediumseagreen", "greenyellow", "gold", "darkorange", "firebrick1")
varit2 <- c("dodgerblue2", "slateblue", "darkorchid", "violetred", "brown1")
varit3 <- c("dodgerblue2", "slateblue", "grey70", "violetred", "brown1")
varit4 <- c("mediumseagreen", "greenyellow", "grey70", "darkgoldenrod1", "darkorange")


#################### Fig 2. Principal coordinates analysis

#Preprocess data
survey <- survey1[survey1$Allhypoanswers , ]
answ <- 6 - sapply(survey1[survey1$Allhypoanswers , allhypos], as.numeric)
colnames(answ) <- as.character(survey1_questions$Brief[survey1_questions$Name %in% allhypos])

pcoa_caps <- capscale(t(answ) ~ 1, distance="euclidean") ##PCoA done
pcoa_caps_all <- pcoa_caps

# traits could be derived from survey1_questions, too.
traits <- as.factor(c(rep("bipedalism", 10), rep("brain", 10), rep("hairlessness", 8),
                      rep("fat", 4), rep("larynx", 3), rep("speech", 7), rep("other", 9)))

## All experts and all hypotheses, data from individuals that answered all hypotheses in one figure
if(FALSE){
colstr <- c("palevioletred1","royalblue1","seagreen1","violet","khaki2","skyblue", "orange")
trait.cols <- colstr[traits]

hypo_sizes <- colMeans(answ)
leg_sizes <- c(4, 3, 2, 1, 0.01)

pdf(file="pcoa_plot.pdf", height=6, width=7.5)
plot(pcoa_caps, display = c("sp", "wa"), type="n", xlim=c(-6,4.5)) ## PCoA biplot, full scale
points(pcoa_caps, display= c("sp"), col="gray40") # adding the people points
points(pcoa_caps, display= c("wa"), pch=19, cex=hypo_sizes, col=trait.cols)
text(pcoa_caps, display=c("wa"), srt=25, cex=0.5)
legend(x=-6, y=3.8, levels(traits), fill=colstr, bty="n", cex=1)
legend(-6, -2, legend=c("Very likely", "Moderately likely",
                        "No opinion", "Moderately unlikely", "Very unlikely"),
       pch=21, pt.cex = leg_sizes, bty="n", cex=1)
dev.off()
}


##### Figure 2: PCoA by expertise

colstr <- c("indianred","royalblue1","seagreen1","violet","khaki2","skyblue", "orange")
trait.cols <- colstr[traits]

hypo_sizes <- colMeans(answ)*0.8
leg_sizes <- c(4, 3.2, 2.4, 1.6, 0.8)

pdf(file="pcoakuvat.pdf", height=10.5, width=7)
par(mar=c(3,3,2,1), mgp=c(1.8,0.7,0))
layout(matrix(c(1,6,2,3,4,5), nrow=3, ncol=2, byrow = TRUE))

xy <- data.frame(scores(pcoa_caps)$sites)
xy2 <- data.frame(scores(pcoa_caps)$species)

# plot(pcoa_caps, display = c("sp", "wa"), type="n", main="All answers") ## PCoA biplot, full scale
plot(xy, xlim=c(-5,4), ylim=c(-5,4), type="n", asp=1, cex.axis=1.2, cex.lab=1.4) ## PCoA biplot, full scale
title(main="All answers", cex.main=1.8)
points(xy2[,1], -xy2[,2], col="gray40", pch=19) # adding the people points
points(xy[,1], -xy[,2], pch=19, cex=hypo_sizes, col=trait.cols)
text(xy[,1], -xy[,2], rownames(xy), srt=25)#, cex=0.8)
legend("topleft", "(a)", cex=1.6, bty="n")

panel <- c("anthropologist", "biologist", "human biologist", "other")
title <- c("Anthropologist", "Biologist", "Human biologist", "Other")

for(i in 1:4){
  
  survey <- survey1[survey1$Allhypoanswers , ]
  survey <- survey[survey$Expertise == panel[i],]
  answ <- 6 - sapply(survey[, allhypos], as.numeric)
  colnames(answ) <- as.character(survey1_questions$Brief[survey1_questions$Name %in% allhypos])
  
  pcoa_caps <- capscale(t(answ) ~ 1, distance="euclidean") ##PCoA done
  
  hypo_sizes <- colMeans(answ)*0.8
  
  xy <- data.frame(scores(pcoa_caps)$sites)
  if(panel[i] =="human biologist") {xy <- -xy}
  if(panel[i] =="other") {xy[,2] <- -xy[,2]}
  xlab <- names(xy)[1]
  ylab <- names(xy)[2]
  
  traits <- as.factor(c(rep("bipedalism", 10), rep("brain", 10), rep("hairlessness", 8),
                        rep("fat", 4), rep("larynx", 3), rep("speech", 7), rep("other", 9)))
  colstr <- c("indianred","royalblue1","seagreen1","violet","khaki2","skyblue", "orange")
  trait.cols <- colstr[traits]
  
  plot(xy[,1], xy[,2], type="p", pch=19, cex=hypo_sizes, col=trait.cols, 
       asp=1, xlab=xlab, ylab=ylab, cex.axis=1.2, cex.lab=1.4)
  title(main=title[i], cex.main=1.8)
  
  text(xy, rownames(xy), srt=25)#, cex=0.6)
  # text(x=c(-3.2,-4,-2.4,-2.7)[i], y=c(2.6,2.7,1.6,1.9)[i], c("b","c","d","e")[i], cex=3)
  legend("topleft", paste("(",letters[i+1],")",sep=""), cex=1.6, bty="n")
  
}

plot(xy, type="n", axes=FALSE, xlab="", ylab="")
legend("top", levels(traits), fill=colstr, bty="n", cex=1.2)
legend("bottom", legend=c("Very likely", "Moderately likely",
                              "No opinion", "Moderately unlikely", "Very unlikely"),
       pch=21, pt.cex = leg_sizes, bty="n", cex=1.2, y.intersp=1.8)

dev.off()

#################### Fig 3. Kuvat ihmisistä
## Something wrong with code

xy <- data.frame(scores(pcoa_caps_all)$species) 

# making colour palette for values 1.0, 1.1, ... , 5.0
pal41 <- colorRampPalette(c("indianred1", "royalblue3"))(41) # 41-step colour slider
# pal41 <- colorRampPalette(c("red2", "royalblue1"))(41) # 41-step colour slider
# pal5 <- c("seagreen1","khaki2","darkorange1","darkblue","black")
# pal3 <- c("seagreen1","seagreen4","darkblue")
pal3 <- c("seagreen1","darkorange1","darkblue")
pal2 <- c("darkblue","seagreen1")
# pal4a <- c("seagreen1","green3","darkblue","purple")
pal4a <- c("seagreen1","khaki2","darkorange1","darkblue")
# pal4b <- c("brown3", "limegreen", "blue", "orange4")
pal4b <- c("blue3", "limegreen", "goldenrod1", "orange4")

pch <- 20
cex <- 0.8
xlim <- c(-0.6,0.5)
ylim <- c(-0.6,0.5)
xlab <- "MDS 1"
ylab <- "MDS 2"
legpos <- c(-0.05,0.55)

pdf(file="peopleplot.pdf", height=8, width=6)
par(mfrow=c(3, 2), mar=c(3,3,2,1), mgp=c(1.6,0.6,0))

# AAH average
plot(xy[,1], xy[,2], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, pch=pch, cex=cex, 
     col=pal41[survey$aah_average*10-9], asp=1,
     main="A) Average AAH score")
legend(legpos[1], legpos[2], levels(survey$Wading), fill=pal41[(5:1)*10-9], bty="n")

# Top12 average
plot(xy[,1], xy[,2], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, pch=pch, cex=cex, 
     col=pal41[survey$nonaah.average12*10-9], asp=1,
     main="B) Average Top12 score")
legend(legpos[1], legpos[2], levels(survey$Wading), fill=pal41[(5:1)*10-9], bty="n")

# Publications
plot(xy[,1], xy[,2], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, pch=pch, cex=cex, 
     col=pal4a[survey$Pub.all.peer], asp=1,
     main="C) Number of publications")
points(aggregate(xy$MDS1, survey["Pub.all.peer"], mean)$x,
       aggregate(xy$MDS2, survey["Pub.all.peer"], mean)$x,
       col="black", pch=3, cex=2, lwd=4)
points(aggregate(xy$MDS1, survey["Pub.all.peer"], mean)$x,
       aggregate(xy$MDS2, survey["Pub.all.peer"], mean)$x,
       col=pal4a, pch=3, cex=1.5)
legend(legpos[1], legpos[2], levels(survey$Pub.all.peer), fill=pal4a, bty="n")

# Expertise
plot(xy[,1], xy[,2], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, pch=pch, cex=cex, 
     col=pal4b[survey$Expertise], asp=1,
     main="D) Field of expertise")
points(aggregate(xy$MDS1, survey["Expertise"], mean)$x,
       aggregate(xy$MDS2, survey["Expertise"], mean)$x,
       col="black", pch=3, cex=2, lwd=4)
points(aggregate(xy$MDS1, survey["Expertise"], mean)$x,
       aggregate(xy$MDS2, survey["Expertise"], mean)$x,
       col=pal4b, pch=3, cex=1.5)
legend(legpos[1], legpos[2], levels(survey$Expertise), fill=pal4b, bty="n")

# Familiarity
plot(xy[,1], xy[,2], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, pch=pch, cex=cex, 
     col=pal3[survey$Familiarity], asp=1,
     main="E) Familiarity with human evolution")
points(aggregate(xy$MDS1, survey["Familiarity"], mean)$x,
       aggregate(xy$MDS2, survey["Familiarity"], mean)$x,
       col="black", pch=3, cex=2, lwd=4)
points(aggregate(xy$MDS1, survey["Familiarity"], mean)$x,
       aggregate(xy$MDS2, survey["Familiarity"], mean)$x,
       col=pal3, pch=3, cex=1.5)
legend(legpos[1], legpos[2], c("Not at all", "Some idea", "Know well"), fill=pal3, bty="n")

# Teaching
plot(xy[,1], xy[,2], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, pch=pch, cex=cex, 
     col=pal2[survey$Teaching], asp=1,
     main="F) Taught human evolution")
points(aggregate(xy$MDS1, survey["Teaching"], mean)$x,
       aggregate(xy$MDS2, survey["Teaching"], mean)$x,
       col="black", pch=3, cex=2, lwd=4)
points(aggregate(xy$MDS1, survey["Teaching"], mean)$x,
       aggregate(xy$MDS2, survey["Teaching"], mean)$x,
       col=pal2, pch=3, cex=1.5)
legend(legpos[1]+0.2, legpos[2], c("Yes", "No"), fill=pal2, bty="n")

dev.off()

############################## Fig 5. Credibility frequences of hypotheses by expertise
### Fig 5 is first because it uses the original ordering of Hypotheses

#darkbrown <- "#291500"
#browngradient <- colorRampPalette(c(darkbrown, "tan3", "lemonchiffon")) # making the gradient for colors
#browns <- browngradient(9)
#browns2 <- browngradient(7)
#browns3 <- browngradient(6)

blues <- rev(c(
  "#f7fcf0",
  "#e0f3db",
  "#ccebc5",
  "#a8ddb5",
  "#7bccc4",
  "#4eb3d3",
  "#2b8cbe",
  "#08589e",
  "black"
))
browns <- rev(c(
  "#ffffe5",
  "#fff7bc",
  "#fee391",
  "#fec44f",
  "#fe9929",
  "#ec7014",
  "#cc4c02",
  "#8c2d04",
  "black"
))

lintyp <- rep(c("solid", "twodash","longdash"), 4)

# Magenta is used for hypotheses where anthropoligists differ from others.
biped.cols <- c(browns[1:7], blues[4], browns[8:9])
brain.cols <- c(browns[1:3], "magenta", browns[5], blues[c(2,4)], browns[7:9])
hair.cols <- c(browns[1], "magenta", browns[3:4], blues[4], browns[6:8])

other <- ggplot(melted_survey[melted_survey$Trait == "Other" , ],			# Other water traits
                aes(x = as.numeric(Answer), colour = Hypothesis))+#, linetype=Hypothesis)) +
  scale_colour_manual(values=blues) + 
  scale_linetype_manual(values = lintyp)+
  geom_freqpoly(binwidth=1, size = 1) +
  coord_cartesian(xlim= c(1.2,4.8)) +
  labs(x = "very likely on left", y = "count", title = "Other") +
  theme(
    legend.position = "bottom", 
    legend.direction = "vertical", 
    axis.text.x = element_blank(),
    legend.key.width = unit(1.2, "cm")
  ) +
  facet_wrap(~ Expertise, scales="free_y", ncol=1)

biped <- ggplot(melted_survey[melted_survey$Trait == "Bipedalism" , ],			#Bipedalism hypotheses
                aes(x = as.numeric(Answer), colour = Hypothesis))+#, linetype=Hypothesis)) +
  scale_colour_manual(values=biped.cols) + 
  scale_linetype_manual(values = lintyp)+
  geom_freqpoly(binwidth=1, size = 1) +
  coord_cartesian(xlim= c(1.2,4.8)) +
  labs(x = "very likely on left", y = "count", title = "Bipedalism") +
  theme(
    legend.position = "bottom", 
    legend.direction = "vertical", 
    axis.text.x = element_blank(),
    legend.key.width = unit(1.2, "cm")
  ) +
  facet_wrap(~ Expertise, scales="free_y", ncol=1)

brain <- ggplot(melted_survey[melted_survey$Trait == "Big brain" , ], 			# Brain hypotheses
                aes(x = as.numeric(Answer), colour = Hypothesis))+#, linetype=Hypothesis)) + 
  scale_colour_manual(values=brain.cols) + 
  scale_linetype_manual(values = lintyp)+
  geom_freqpoly(binwidth=1, size = 1) +
  coord_cartesian(xlim= c(1.2,4.8)) +
  labs(x = "very likely on left", y = "count", title = "Big brain") +
  theme(
    legend.position = "bottom", 
    legend.direction = "vertical", 
    axis.text.x = element_blank(),
    legend.key.width = unit(1.2, "cm")
  ) +
  facet_wrap(~ Expertise, scales="free_y", ncol=1)

hair <- ggplot(melted_survey[melted_survey$Trait == "Nakedness" , ],			# Hairlessness hypotheses
               aes(x = as.numeric(Answer), colour = Hypothesis))+#, linetype=Hypothesis)) + 
  scale_colour_manual(values=hair.cols) + 
  scale_linetype_manual(values = lintyp)+
  geom_freqpoly(binwidth=1, size = 1) +
  coord_cartesian(xlim= c(1.2,4.8)) +
  labs(x = "very likely on left", y = "count", title = "Nakedness") +
  theme(
    legend.position = "bottom", 
    legend.direction = "vertical", 
    axis.text.x = element_blank(),
    legend.key.width = unit(1.2, "cm")
  ) +
  facet_wrap(~ Expertise, scales="free_y", ncol=1)

pdf(file="lineplot.pdf", height=11, width=8)
grid.arrange(biped, brain, hair, other, ncol=4)
dev.off()

############# Fig 4. THE GEOM_BAR TO SHOW THE PROPORTIONS OF ANSWERS FOR ALL HYPOTHESES

## Muuttaa hypoteesien järjestyksen vastakkaiseksi. Muita kuvia tehdessä halutaan, että ne on 
# päinvastaisessa järjestyksessä, joten pitää huomata tämä muita piirrettäiessä.
melted_survey$Hypothesis <- factor(
  melted_survey$Hypothesis,
  levels = rev(levels(melted_survey$Hypothesis)), ordered=TRUE)


pdf(file="all_answers_green.pdf", height=6, width=7.5)
ggplot(melted_survey, aes(Hypothesis, fill=Answer)) +		## Plotted by trait
  geom_bar(position="fill") +
  scale_fill_manual(values=varit4) +
  labs(x="", y="Percentage of respondents", title="Support for different hypotheses") +
  theme_grey(base_size=9) +
  theme(axis.text.x=element_blank(), plot.title = element_text(hjust = 0),
        legend.position = "bottom", legend.direction = "horizontal") + #  legend.text.align = 0,
  coord_flip() +
  scale_y_reverse() + 
  facet_wrap(~Trait, ncol=2, scales="free")
dev.off()

#lisäämällä tämän scale_fill_manualin sisälle saa legendistä titlen pois ja asetettua ne kahteen sarakkeeseen
# guide = guide_legend(title=NULL, label.position = "bottom", nrow=3, byrow=TRUE,
# title.hjust = 1.3)

########################### Fig 6. Number of credible hypotheses

hypos <- allhypos[-(71:79-28)] #"Other" trait removed because not meaningful here.
traits <- as.character(survey1_questions$Trait[
  match(
    colnames(survey1[hypos]), 
    survey1_questions$Name)
  ])

numHypotheses <- function(credibility) {
  # support is an array for tests whether a trait*respondent pair
  # fulfils the given credibility criteria (1:very likely, 5: very unlikely).
  support <- sapply(survey1[hypos], FUN = as.numeric) <= 6-credibility
  colnames(support) <- gsub("\\.", " ", colnames(support))
  
  # test is an array with number of credible hypothesis for each 
  # respondent * trait pair.
  test <- sapply(
    unique(traits), FUN = function(x) {
      rowSums(support[,traits %in% x], na.rm = TRUE)
    }
  )
  # Combine with the determinant columns in survey1 and melt.
  test <- cbind(survey1[1:28], test)
  test <- melt(
    test, 
    measure.vars = unique(traits), 
    variable.name = "Trait", 
    value.name = "Credible.hypotheses"
  )
  return(test)
}

numhypos <- rbind(
  data.frame(numHypotheses(5), Likelihood = "Very likely"),
  data.frame(numHypotheses(4), Likelihood = "Likely")
)
numhypos$Credible.hypotheses <- factor(
  numhypos$Credible.hypotheses, 
  levels = as.character(0:10),
  ordered = TRUE
)

numhypos <- merge(
  numhypos, 
  aggregate(survey1_questions$Name, by = survey1_questions["Trait"], length)
)
numhypos$Trait2 <- paste(numhypos$Trait," (",numhypos$x,")", sep="")
numhypos$Trait2 <- factor(numhypos$Trait2, levels=rev(sort(unique(numhypos$Trait2))))

# bipedalism (10), encephalization (10), nakedness (8), descended larynx (3),
# subcutaneous fat (4), and speech (7).

pdf(file="likely_hypotheses.pdf", height=5, width=7)
ggplot(numhypos, aes(x=Trait2, fill=Credible.hypotheses))+
  geom_bar(position="fill")+coord_flip()+ facet_wrap(~ Likelihood)+
  scale_fill_brewer(palette = "Set3", name="Number of\nhypotheses")+
  labs(title = "Counting likely or very likely hypotheses",
       y = "Percentage of respondents", x = "Trait")+
  theme(axis.text.y = element_text(size=10), 
        legend.text= element_text(size=10),
        axis.text.x = element_blank())+
  scale_y_reverse()
dev.off()

############## Fig 7. AAH criticism in the expertise groups

pdf(file="aah_criticism_green.pdf", height=5, width=7.5)
ggplot(critique, aes(x = Critique, fill=Answer)) +	
  geom_bar(position="fill") +
  scale_fill_manual(values=varit4) + facet_grid(. ~ Expertise)+
  labs(x="", y="Percentage of respondents", title="Support for AAH criticism") +
  theme_grey(base_size=9) +
  theme(axis.text.x=element_blank(), #plot.title = element_text(hjust = 1.3), 
        legend.position = "bottom", legend.text.align = 0, legend.direction = "horizontal") +
  coord_flip() +
  scale_y_reverse()
dev.off()

############# Supplementary figures: credibility score by publication numbers.
## Something wrong with code

temp <- melted_survey
temp$Answer <- 6 - as.numeric(temp$Answer)
nam <- gsub("\\.", " ", colnames(survey1)[29:79])
temp$Hypothesis <- factor(temp$Hypothesis, levels = nam)

temp2 <- na.omit(temp[c("Answer","Pub.all.peer","Hypothesis")])
colnames(temp2)[colnames(temp2)=="Answer"]<- "Result"
levels(temp2$Pub.all.peer)[levels(temp2$Pub.all.peer) %in% c("none", "1-10")] <- "0-10"
temp2 <- Ovariable(output=temp2, marginal=colnames(temp2) %in% c("Pub.all.peer","Hypothesis"))
temp2 <- summary(temp2, fun=c("mean", "sd", "length"))
temp2$ci <- temp2$sd/sqrt(temp2$length)*1.96

pdf(file="credibility score by all publications.pdf", height=12, width=17.5)
ggplot(temp2, aes(x = Pub.all.peer, y = mean, group=Hypothesis))+facet_wrap(~Hypothesis)+
  geom_errorbar(aes(ymin=mean-ci, ymax=mean+ci, width = 0.1))+
  geom_line()+
  geom_point()+
  labs(x= "Number of all scientific publications", y = "Credibility score +- 95% CI")
dev.off()

temp2 <- na.omit(temp[c("Answer","Pub.hum.peer","Hypothesis")])
colnames(temp2)[colnames(temp2)=="Answer"]<- "Result"
temp2 <- Ovariable(output=temp2, marginal=colnames(temp2) %in% c("Pub.hum.peer","Hypothesis"))
temp2 <- summary(temp2, fun=c("mean", "sd", "length"))
temp2$ci <- temp2$sd/sqrt(temp2$length)*1.96

pdf(file="credibility score by hum_evol publications.pdf", height=12, width=17.5)
ggplot(temp2, aes(x = Pub.hum.peer, y = mean, group=Hypothesis))+facet_wrap(~Hypothesis)+
  geom_errorbar(aes(ymin=mean-ci, ymax=mean+ci, width = 0.1))+
  geom_line()+
  geom_point()+
  labs(x= "Number of scientific publications about human evolution", y = "Credibility score +- 95% CI")
dev.off()

How scientists perceive the evolutionary origin of human traits: results of a survey study

Contents

Article

Materials and methods

Data

Calculations

Preprocess for Dryad data

Correlation and histograms

Principal coordinates and component analyses, linear regression

AAH criticism and anova

Images

See also

Navigation menu