Last active
April 11, 2016 07:53
-
-
Save MartinMacharia/a233f895da55a023a57a to your computer and use it in GitHub Desktop.
Mozambique OFRA Workshop, quick mention about R
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
yield=read.csv(file.choose(),header=T) | |
yield | |
#Another way of importing the data | |
#SmallYieldData <- read.csv("C:/Users/Administrator/Desktop/SmallYieldData.csv") | |
#View(SmallYieldData) | |
#You can also type the data in R as follows | |
yield=data.frame(Nrates=c(0,10,20,30,40,50,60),yield=c(2748,3162,3307,3571,3576,3544,3600)) | |
yield | |
#You can view in which directory you are in R and also change the working direcory | |
#as follows | |
getwd() | |
setwd("C:/Users/Administrator/Desktop") | |
SmallYieldData | |
#Kandara, 2014, Bean Season A | |
Kandara_Bean=read.csv(file.choose(), header=T) | |
Kandara_Bean | |
summary(Kandara_Bean) | |
attach(Kandara_Bean) | |
#Vary P and Krates for different plots | |
par(mfrow=c(1,2)) | |
boxplot(Yield~Nrate, subset=Prate==22.5 & Krate==0) | |
plot(Yield~Nrate, subset=Prate==22.5 & Krate==0) | |
par(mfrow=c(1,2)) | |
boxplot(Yield~Nrate, subset=Prate==0 & Krate==0 & Manure==0) | |
plot(Yield~Nrate, subset=Prate==0 & Krate==0 & Manure==0) | |
Kandara2014bean=read.csv(file.choose(), header=T) | |
Kandara2014bean | |
Assypmodel=nls(Mean~a-b*(c^Nrate), start=list(a=2000,b=1000, c=0.9), data=Kandara2014bean) | |
Assypmodel | |
summary(Assypmodel) | |
#Anova | |
summary(Kandara_Bean) | |
str(Kandara_Bean) | |
Nrate1=as.factor(Nrate) | |
Prate1=as.factor(Prate) | |
Krate1=as.factor(Krate) | |
Manure1=as.factor(Manure) | |
Diagnostic1=as.factor(Diagnostic) | |
Rep1=as.factor(Rep) | |
kandarabean1=data.frame(Rep1,Nrate1,Prate1,Krate1,Manure1,Diagnostic1,Yield) | |
kandarabean1 | |
str(kandarabean1) | |
fit=aov(Yield~Rep1+Nrate1+Prate1+Krate1+Manure1+Diagnostic1, data=kandarabean1) | |
fit | |
summary(fit) | |
fit1=aov(Yield~Rep1+Nrate1+Prate1+Krate1, data=kandarabean1) | |
fit1 | |
summary(fit1) | |
library(lme4) | |
fit2=lmer(Yield~Nrate1+Prate1+Krate1+(1|Rep1), data=kandarabean1) | |
fit2 | |
summary(fit2) | |
#As a single treatment effect | |
kandarabean2=read.csv(file.choose(), header=T) | |
kandarabean2 | |
str(kandarabean2) | |
summary(kandarabean2) | |
attach(kandarabean2) | |
Trt1=as.factor(Trt) | |
Rep1=as.factor(Rep) | |
kandarabean3=data.frame(Rep1,Trt1,Yield) | |
kandarabean3 | |
boxplot(Yield~Trt,data=kandarabean3) | |
fit3=aov(Yield~Rep1+Trt1, data=kandarabean3) | |
fit3 | |
summary(fit3) | |
plot(fit3) | |
TukeyHSD(fit3) | |
print(model.tables(fit3,"means"),digits=3) | |
#KandaraMaize2014 | |
Kandaramaize2014=read.csv(file.choose(), header=T) | |
Kandaramaize2014 | |
attach(Kandaramaize2014) | |
boxplot(Yield~Trt)#, subset=Prate==0 & Krate==0 & Manure==0 & Diagnosis==0) | |
plot(Yield~Nrate, subset=Prate==0 & Krate==0 & Manure==0 & Diagnosis==0) | |
newdata2 <- subset(Kandaramaize2014, Prate==0 & Krate==0 & Manure==0 & Diagnosis==0, | |
select=c(Yield, Nrate,Prate,Krate,Manure,Diagnosis)) | |
newdata2 | |
str(newdata2) | |
plot(Yield~Nrate, data=newdata2) | |
Mod1=nls(Mean~a-b*(c^Nrate), start=list(a=3500,b=500,c=0.5), data=newdata2) | |
Mod | |
################################################### | |
#Kisii bean 2014 | |
kisiibean2014=read.csv(file.choose(),header=T) | |
kisiibean2014 | |
names(kisiibean2014) | |
plot(Yield~Nrate, subset=Prate==0 & Krate==0 & Manure==0 & Diagnosis==0, data=kisiibean2014) | |
newdata3 <- subset(kisiibean2014, Prate==0 & Krate==0 & Manure==0 & Diagnosis==0, | |
select=c(Yield, Nrate,Prate,Krate,Manure,Diagnosis)) | |
newdata3 | |
Mod3=nls(Yield~a-b*(c^Nrate), start=list(a=3500,b=2500,c=0.9), data=newdata3) | |
Mod3 | |
kisii2014beanN=read.csv(file.choose(),header=T) | |
kisii2014beanN | |
plot(Yield~Nrate, data=kisii2014beanN) | |
Mod4=nls(Yield~a-b*(c^Nrate), start=list(a=0.9,b=0.3,c=0.9), data=kisii2014beanN) | |
Mod4 | |
############################################################################## | |
N=c(0,10,20,30) | |
yield=c(0.550925926,0.760277778,0.768425926,0.768425926) | |
dat <- as.data.frame(cbind(N,yield)) | |
dat | |
plot(yield~N, data=dat) | |
mod9=nls(yield~a-b*(c^N), start=list(a=0.3, b=0.8,c=0.9),data=dat) | |
mod9 | |
summary(mod9) | |
coef(mod9) | |
#############linear to plateau | |
f.lrp=function(N,a,b,t.x){ | |
ifelse(N>t.x,a+b*t.x,a+b*N) | |
} | |
#t.x is maximum fertilizer that will give any result | |
m <- nls(yield ~ f.lrp(N, a, b, t.x), data =dat, start = list(a = 0.5, b = 0.2, t.x = 10), trace = T ) | |
m | |
summary(m) | |
coefficients(m) | |
############################################################# | |
#Covariate analysis | |
CovKisii2014=read.csv(file.choose(),header=T) | |
CovKisii2014 | |
str(CovKisii2014) | |
names(CovKisii2014) | |
Rep1=as.factor(CovKisii2014$Rep) | |
str(Rep1) | |
Rep1 | |
CovKisii2014$Rep | |
Nrate1=as.factor(CovKisii2014$Nrate) | |
Nrate1 | |
str(Nrate1) | |
Yield1=CovKisii2014$Yield | |
Count1=CovKisii2014$Plant.count | |
newcov=data.frame(Rep1,Nrate1,Yield1,Count1) | |
newcov | |
str(newcov) | |
#mixed model approach | |
library(lme4) | |
fit3=lmer(Yield1~Nrate1+Count1+(1|Rep1), data=newcov) | |
fit3 | |
summary(fit3) | |
#anova appoach-Analysis of Covariance (ANCOVA) | |
#To test for differences between different mean Nrates when | |
#we know that an extraneous continuous variable (Plant count) | |
#affects the continuous outcome variable (Yield) | |
fit5=aov(Yield1~Nrate1+Rep1+Count1,data=newcov) | |
fit5 | |
summary(fit5) | |
print(model.tables(fit5 ,"means"),digits=3) | |
#However, note that aov() will use Type I Sum of Squares by | |
#default. So our result here will depend upon the order in | |
#which we entered the predictors. | |
#Better to use Type III Sum of Squares, available using the | |
#Anova() function of the car package: | |
require(car) | |
###################################################################### | |
fit6=aov(Count1~Nrate1+Rep1, data=newcov) | |
summary(fit6) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
This was a quick introduction to R, future followup work shall be posted here