Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save MartinMacharia/a233f895da55a023a57a to your computer and use it in GitHub Desktop.
Save MartinMacharia/a233f895da55a023a57a to your computer and use it in GitHub Desktop.
Mozambique OFRA Workshop, quick mention about R
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)
@MartinMacharia
Copy link
Author

This was a quick introduction to R, future followup work shall be posted here

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment