1 / 20

Exercise 3 week 2

Exercise 3 week 2. Var 1. Var 2. Experimental setup. Climate 2. Climate 1. Fert 1. Fert 2. Fert 3. Block 1. Var 1. Var 2. Var 1. Var 2. Block 2. Block 3. Section. Variables. Response variable: End level Systematic part: Climate, Fertilizer, Variaty Random effects Block,

nona
Télécharger la présentation

Exercise 3 week 2

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Exercise 3 week 2

  2. Var 1 Var 2 Experimental setup Climate 2 Climate 1 Fert 1 Fert 2 Fert 3 Block 1 Var 1 Var 2 Var 1 Var 2 Block 2 Block 3 Section

  3. Variables • Response variable: • End level • Systematic part: • Climate, Fertilizer, Variaty • Random effects • Block, • Section, • Plot (factor(section):factor(fertilizer)

  4. Endlevel= a(climatei, fertilizeri, varietyi)+ d(blocki),f(sectioni)+ g(ploti) + ei • i = 1…. 36 • ei ~ N(0,d2) • d ~N(0, dd2), e ~N(0, de2), g ~N(0, dg2) • plot = factor(section):factor(fertilizer)

  5. Reading in data data=read.table(file.choose(), header=T) library(nlme) data$Plot<-with(data,factor(section):factor(fert)) _____________________________________________________________ > head(data,3) block section clim fert variety endlevel rate Plot 1 1 1 2 2.0 aminex 48.8981 0.06915 1:2 2 1 1 2 2.0 dalibor 42.2463 0.06595 1:2 3 1 1 2 3.5 aminex 48.2108 0.04679 1:3.5 _________________________________________________________ is.factor(data$block) # Test to see if the variables are factors or covariates is.factor(data$Plot) is.factor(data$clim) is.factor(data$variety) is.factor(data$fert) Is.factor(section) data$Blok<-with(data,factor(block)) # Define block as a factor which can be used in lme data$Section<-with(data,factor(section)) # Define section as a factor which can be used in lme

  6. Start model model1= lme(endlevel ~ factor(clim)*factor(fert)*factor(variety), random=~1|Blok/Section/Plot, method="REML", data=data) ############### model validation lm on start model: model1.ml model1ml= lm(endlevel ~ factor(clim)*factor(fert)*variety+Plot, data=data) par(mfrow=c(1,2)) pred<-fitted(model1ml) raw.res<-residuals(model1ml) plot(pred,raw.res) # Predicted values vs. Raw residuals abline(h=c(-1.96*sd(raw.res), 0, 1.96*sd(raw.res))) qqnorm(raw.res) # Standardized residuals vs. Quintiles # of standardized normal distribution qqline(raw.res)

  7. Model validation on model: model1 fitted with lm # This model validation looks quite ok

  8. Model validation lm model1ml= lm(endlevel ~ factor(clim)*factor(fert)*variety*Blok*Plot*Section, data=data) pred<-fitted(model1lm) raw.res<-residuals(model1lm) plot(pred,raw.res) abline(h=c(-1.96*sd(raw.res), 0, 1.96*sd(raw.res))) Make instead model1ml= lm(endlevel ~ factor(clim) *factor(fert)*variety+Plot Blok and Section is nested in Plot

  9. Question Why does the EBLUPS give a qqnrom plot with a plateau ? model1= lme(endlevel ~ factor(clim)*factor(fert)*factor(variety), random=~1|Blok/Section/Plot, method="REML", data=data) qqnorm(model1) eblups<-as.vector(unlist(ranef(model1))) qqnorm(eblups) abline(0,sd(eblups))

  10. Start model validation EBLUPS par(mfrow=c(1,3)) eblups.Blok<-as.vector(unlist(ranef(model1c)["Blok"])) qqnorm(eblups.Blok); abline(0,sd(eblups.Blok)) eblups.Section<-as.vector(unlist(ranef(model1c)["Section"])) qqnorm(eblups.Section);abline(0,sd(eblups.Section)) eblups.Plot<-as.vector(unlist(ranef(model1c)["Plot"])); qqnorm(eblups.Plot);abline(0,sd(eblups.Plot))

  11. Summary model 1 Linear mixed-effects model fit by REML Data: data AIC BIC logLik 188.3954 207.2442 -78.19769 Random effects: Formula: ~1 | Blok (Intercept) StdDev: 3.387274 Formula: ~1 | Section %in% Blok (Intercept) StdDev: 2.074986 Formula: ~1 | Plot %in% Section %in% Blok (Intercept) Residual StdDev: 0.0008466898 4.175731 Fixed effects: endlevel ~ factor(clim) * factor(fert) * variety Value Std.Error DF t-value p-value (Intercept) 50.95107 3.327458 12 15.312310 0.0000 factor(clim)2 0.73287 3.807213 2 0.192494 0.8651 factor(fert)3.5 -3.28387 3.409470 8 -0.963160 0.3637 factor(fert)4 4.13687 3.409470 8 1.213346 0.2596 varietydalibor -11.98480 3.409470 12 -3.515150 0.0043 factor(clim)2:factor(fert)3.5 -1.75333 4.821719 8 -0.363632 0.7256 factor(clim)2:factor(fert)4 -5.45873 4.821719 8 -1.132113 0.2904 factor(clim)2:varietydalibor 2.03050 4.821719 12 0.421115 0.6811 factor(fert)3.5:varietydalibor 5.35603 4.821719 12 1.110814 0.2884 factor(fert)4:varietydalibor -5.74887 4.821719 12 -1.192286 0.2562 factor(clim)2:factor(fert)3.5:varietydalibor 2.09800 6.818941 12 0.307672 0.7636 factor(clim)2:factor(fert)4:varietydalibor 6.69730 6.818941 12 0.982161 0.3454

  12. Model reduction climat=factor(data$clim);ferti=factor(data$fert) model1ML= lme(endlevel ~ climat*ferti*variety , random=~1|Blok/Section/Plot, method="ML", data=data) # Start model model2ML= lme(endlevel ~ climat*ferti+climat*variety+ferti*variety , random=~1|Blok/Section/Plot, method="ML", data=data) # model2ML taken out Triple effect modification model3ML= lme(endlevel ~ climat+ferti + climat*variety + ferti*variety , random=~1|Blok/Section/Plot, method="ML", data=data) # model3ML taken out effect modification clim : fert ____________________________________________________________________ anova(model1ML,model2ML)# P-value 0.4778 NS # Triple effect modification Possible anova(model2ML,model3ML)# P-value 0.7551 NS # clim : fert Possible ____________________________________________________________________

  13. Model reduction continued model4ML= lme(endlevel ~ climat + ferti + variety + climat: variety , random=~1|Blok/Section/Plot, method="ML", data=data) # model4ML Taken out effect modification fert : variety model5ML= lme(endlevel ~ climat + ferti + variety + ferti:variety , random=~1|Blok/Section/Plot, method="ML", data=data) # model5ML taken out effect modification , clim : variety ___________________________________________________ anova(model3ML,model4ML)# P-value 0.0134 * # clim : fert , fert : variety Not Possible anova(model3ML,model5ML)# P-value 0.0416 * # clim : fert , clim : variety Not Possible ___________________________________________________ End model is model: model3ML

  14. According to summary end model model3ML reduction is still possible model3REML= lme(endlevel ~ climat+ferti+climat*variety+ferti*variety , random=~1|Blok/Section/Plot, method="REML", data=data) anova(model3REML) # numDF denDF F-value p-value #(Intercept) 1 14 413.3631 <.0001 #climat 1 2 0.1365 0.7473 #ferti 2 10 0.2982 0.7485 #variety 1 14 53.7166 <.0001 #climat:variety 1 14 3.5624 0.0800 #ferti:variety 2 14 3.9967 0.0423 Deside to focus on climate and take out effect modification climat:variety and climate

  15. Model reduction round 2 model5ML= lme(endlevel ~ climat + ferti + variety + ferti:variety , random=~1|Blok/Section/Plot, method="ML", data=data) # Take out effect modification climate : variaty model6ML= lme(endlevel ~ climat + ferti + variety , random=~1|Blok/Section/Plot, method="ML", data=data) # Take out effect modification ferti:variety model7ML= lme(endlevel ~ ferti * variety , random=~1|Blok/Section/Plot, method="ML", data=data) # Take out systematic effect climate model8ML= lme(endlevel ~ ferti + variety , random=~1|Blok/Section/Plot, method="ML", data=data) # Take out effect modification ferti:variety anova(model5ML,model6ML)#P-value : 0.0219 * Not Possible anova(model5ML,model7ML)#P-value : 0.6563 NS Possible anova(model5ML,model8ML)#P-value : 0.0495 * Not possible

  16. Validation of end model model7lm

  17. Estimates final model summary(model7REML) _________________________________________________________________________________ Fixed effects: endlevel ~ ferti * variety Value Std.Error DF t-value p-value (Intercept) 51.31750 2.721758 15 18.854541 0.0000 ferti3.5 -4.16053 2.390720 10 -1.740284 0.1124 ferti4 1.40750 2.390720 10 0.588735 0.5691 varietydalibor -10.96955 2.390720 15 -4.588387 0.0004 ferti3.5:varietydalibor 6.40503 3.380989 15 1.894426 0.0776 ferti4:varietydalibor -2.40022 3.380989 15 -0.709916 0.4886 ___________________________________________________________________________ model7REML1= lme(endlevel ~ ferti * variety -1, random=~1|Blok/Section/Plot, method="REML", data=data) intervals(model7REML1) _______________________________________________________________________________________ Fixed effects: lower est. upper ferti2 45.2530450 51.317500 57.381955 ferti3.5 41.0925117 47.156967 53.221422 ferti4 46.6605450 52.725000 58.789455 varietydalibor -16.0376506 -10.969550 -5.901449 ferti3.5:varietydalibor -0.7623433 6.405033 13.572410 ferti4:varietydalibor -9.5675933 -2.400217 4.767160 _______________________________________________________________________________________ Question: Can systematic effect -1 be used to report estimates and 95% confidence Int.

  18. Reporting random effects Blok StdDev: 3.537 Section %in% Blok StdDev: 1.509 Plot %in% Section %in% Blok StdDev: 0.000292 Residual StdDev: 4.141 Startmodel Blok StdDev: 3.387 Section %in% Blok StdDev: 2.075 Plot %in% Section %in% Blok StdDev: 0.000847 Residual StdDev: 4.176

  19. Reporting Statistical methods: The combined effects of climate conditions, fertilizer levels and cucumber variety on end level of percent infection was modeled by a linear mixed model with block, section and plot as random effect. Plot was calculated as the effect modification between section and fertilizer treating the variables as factors. Reported p-values correspond to likely hood ratio and were Evaluated at 5% significance level. All analyses were made in R version 2.14.0 (www.r-project.org). Results: There was no triple effect modification of climate conditions, fertilizer levels on end level of percent infection breed (p= 0.48). There was no effect modification of climate conditions and fertilizer level on end level of percent infection (p=0.75). Likewise there seemed to be no effect modification of climate conditions and variety on end level of percent infection (p= 0.080). Climate conditions had no influence on the end level percent infection (p= 0.66). The effect modification between fertilizer level and variety was significant (p= 0.049). Fertilizer level and variety influenced the end level of percent infection. The difference between the two varieties was significant (p-value 0.004) giving an estimate and 95% confidence interval for the difference of -10.97(-16.038; -5.901) for the treatment with fertilizer level 2.0. Furthermore there was no significant differences between the fertilizer levels 3.5 and 4.0 for the variety aminex when compared to the fertilizer level 2.0. This was also the case for the variety dalibor.

More Related