Re: [R] How to extract partial predictions, package mgcv
Thanks Simon, I wasn't aware that predict.gam could be used in this situation. I followed you advice and managed to extract the data I needed. For people interested, I add the code with comments I used here: # # Full code for extracting partial predictions # Start the package mgcv library(mgcv) # Simulate some data (this is from Simon Wood's example in the # package mgcv manual) n-200 sig - 2 dat - gamSim(1,n=n,scale=sig) # Check the data dat ## It looks alright :-) # Run the GAM b-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat) # Plot the partial predictions (You may need to rescale your window to # see the non-linearity par(mfrow=c(1,3)) plot(b, scale=0) ### Clearly, the variables x0 and x2 are non-linear! # I would like to extract the solid line in the graphs and the # associated standard errors from the plots. Thus, I use predict # and change to a data.frame: c-data.frame(predict(b,type=terms,se.fit=TRUE)$fit) names(c)-c(pp_s.x0.,pp_s.I.x1.2..,pp_s.x2.) d-data.frame(predict(b,type=terms,se.fit=TRUE)$se.fit) names(d)-c(se_s.x0.,se_s.I.x1.2..,se_s.x2.) # Merge the three datasets: f=cbind(dat,c,d) #Restrict the number of variables to the ones I am interested in: f-f[,c(x0,pp_s.x0., se_s.x0.)] names(f) ### This data, i.e. f, can now be exported or used for further calculations within R #plot the data par(mfrow=c(1,1)) plot(pp_s.x0.~x0,type=p,pch=2,cex=0.50) sequence=order(x0) lines(x0[sequence],pp_s.x0.[sequence]) lines rug(jitter(x0)) ## Thanks, Staffan -- Staffan Roos, PhD Research Ecologist BTO Scotland Staffan, Take a look at ?predict.gam. You need to use type=terms with se=TRUE to get the quantities that you need. Simon ps. with binary data, method=REML or method=ML for the gam fit are often somewhat more reliable than the default. On Monday 14 September 2009 19:30, Staffan Roos wrote: Dear package mgcv users, I am using package mgcv to describe presence of a migratory bird species as a function of several variables, including year, day number (i.e. day-of-the-year), duration of survey, latitude and longitude. Thus, the global model is: global_model-gam(present ~ as.factor(year) + s(dayno, k=5) + s(duration, k=5) + s(x, k=5) + s(y, k=5), family = binomial(link = logit), data = presence, gamma =1.4) The model works fine, suggesting that all the variables have strong, non-linear, influences on the probability of detecting the species. My main interest is in the effect dayno has on presence, given the inclusion of the other explanatory variables. Thus, I would like to extract the values of the partial prediction of dayno and its associated 2 standard errors above and below the main effect, as shown by the code plot(global_model). I have tried to extract the values by using fitted.values(global_model,dayno), but when plotted, the figure doesn't look like the partial prediction for dayno. Instead, it gives me a very scattered figure (shotgun effect). Has anyone tried to extract the partial predictions? If so, please could you advise me how to do this? Thanks, Staffan P.S.. For comparison, please have a look at Simon Woods paper in R News, 1(2):20-25, June 2001, especially the figures showing the partial predictions of Mackerel egg densities. Using those graphs as an example, I would like to extract the partial predictions for e.g. s(b.depth), given the inclusion of the other variables. -- Staffan Roos, PhD Research Ecologist BTO Scotland __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK +44 1225 386603 www.maths.bath.ac.uk/~sw283 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- View this message in context: http://www.nabble.com/How-to-extract-partial-predictions%2C-package-mgcv-tp25441149p25476107.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to extract partial predictions, package mgcv
It is nice to see worked examples, but there were some errors that relate to not including dataframe references. I've paste in code that runs without error on my machine: On Sep 16, 2009, at 12:39 PM, Staffan Roos wrote: Thanks Simon, I wasn't aware that predict.gam could be used in this situation. I followed you advice and managed to extract the data I needed. For people interested, I add the code with comments I used here: # # Full code for extracting partial predictions # Start the package mgcv library(mgcv) # Simulate some data (this is from Simon Wood's example in the # package mgcv manual) n-200 sig - 2 dat - gamSim(1,n=n,scale=sig) # Check the data dat ## It looks alright :-) # Run the GAM b-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat) # Plot the partial predictions (You may need to rescale your window to # see the non-linearity par(mfrow=c(1,3)) plot(b, scale=0) ### Clearly, the variables x0 and x2 are non-linear! # I would like to extract the solid line in the graphs and the # associated standard errors from the plots. Thus, I use predict # and change to a data.frame: c-data.frame(predict(b,type=terms,se.fit=TRUE)$fit) names(c)-c(pp_s.x0.,pp_s.I.x1.2..,pp_s.x2.) d-data.frame(predict(b,type=terms,se.fit=TRUE)$se.fit) names(d)-c(se_s.x0.,se_s.I.x1.2..,se_s.x2.) # Merge the three datasets: f=cbind(dat,c,d) #Restrict the number of variables to the ones I am interested in: f-f[,c(x0,pp_s.x0., se_s.x0.)] names(f) ### This data, i.e. f, can now be exported or used for further ### calculations within R #plot the data par(mfrow=c(1,1)) plot(pp_s.x0.~x0,type=p,pch=2,cex=0.50, data=f) sequence=order(f$x0) lines(f$x0[sequence],f$pp_s.x0.[sequence]) lines rug(jitter(f$x0)) ## Staffan, Take a look at ?predict.gam. You need to use type=terms with se=TRUE to get the quantities that you need. Simon ps. with binary data, method=REML or method=ML for the gam fit are often somewhat more reliable than the default. On Monday 14 September 2009 19:30, Staffan Roos wrote: Dear package mgcv users, I am using package mgcv to describe presence of a migratory bird species as a function of several variables, including year, day number (i.e. day-of-the-year), duration of survey, latitude and longitude. Thus, the global model is: global_model-gam(present ~ as.factor(year) + s(dayno, k=5) + s(duration, k=5) + s(x, k=5) + s(y, k=5), family = binomial(link = logit), data = presence, gamma =1.4) The model works fine, suggesting that all the variables have strong, non-linear, influences on the probability of detecting the species. My main interest is in the effect dayno has on presence, given the inclusion of the other explanatory variables. Thus, I would like to extract the values of the partial prediction of dayno and its associated 2 standard errors above and below the main effect, as shown by the code plot(global_model). I have tried to extract the values by using fitted.values(global_model,dayno), but when plotted, the figure doesn't look like the partial prediction for dayno. Instead, it gives me a very scattered figure (shotgun effect). Has anyone tried to extract the partial predictions? If so, please could you advise me how to do this? Staffan P.S.. For comparison, please have a look at Simon Woods paper in R News, 1(2):20-25, June 2001, especially the figures showing the partial predictions of Mackerel egg densities. Using those graphs as an example, I would like to extract the partial predictions for e.g. s(b.depth), given the inclusion of the other variables. David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to extract partial predictions, package mgcv
But removing the extraneous second to last line that says just: lines ... would leave the console looking less puzzling. -- David On Sep 16, 2009, at 1:02 PM, David Winsemius wrote: It is nice to see worked examples, but there were some errors that relate to not including dataframe references. I've paste in code that runs without error on my machine: On Sep 16, 2009, at 12:39 PM, Staffan Roos wrote: Thanks Simon, I wasn't aware that predict.gam could be used in this situation. I followed you advice and managed to extract the data I needed. For people interested, I add the code with comments I used here: # # Full code for extracting partial predictions # Start the package mgcv library(mgcv) # Simulate some data (this is from Simon Wood's example in the # package mgcv manual) n-200 sig - 2 dat - gamSim(1,n=n,scale=sig) # Check the data dat ## It looks alright :-) # Run the GAM b-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat) # Plot the partial predictions (You may need to rescale your window to # see the non-linearity par(mfrow=c(1,3)) plot(b, scale=0) ### Clearly, the variables x0 and x2 are non-linear! # I would like to extract the solid line in the graphs and the # associated standard errors from the plots. Thus, I use predict # and change to a data.frame: c-data.frame(predict(b,type=terms,se.fit=TRUE)$fit) names(c)-c(pp_s.x0.,pp_s.I.x1.2..,pp_s.x2.) d-data.frame(predict(b,type=terms,se.fit=TRUE)$se.fit) names(d)-c(se_s.x0.,se_s.I.x1.2..,se_s.x2.) # Merge the three datasets: f=cbind(dat,c,d) #Restrict the number of variables to the ones I am interested in: f-f[,c(x0,pp_s.x0., se_s.x0.)] names(f) ### This data, i.e. f, can now be exported or used for further ### calculations within R #plot the data par(mfrow=c(1,1)) plot(pp_s.x0.~x0,type=p,pch=2,cex=0.50, data=f) sequence=order(f$x0) lines(f$x0[sequence],f$pp_s.x0.[sequence]) rug(jitter(f$x0)) ## Staffan, Take a look at ?predict.gam. You need to use type=terms with se=TRUE to get the quantities that you need. Simon ps. with binary data, method=REML or method=ML for the gam fit are often somewhat more reliable than the default. On Monday 14 September 2009 19:30, Staffan Roos wrote: Dear package mgcv users, I am using package mgcv to describe presence of a migratory bird species as a function of several variables, including year, day number (i.e. day-of-the-year), duration of survey, latitude and longitude. Thus, the global model is: global_model-gam(present ~ as.factor(year) + s(dayno, k=5) + s(duration, k=5) + s(x, k=5) + s(y, k=5), family = binomial(link = logit), data = presence, gamma =1.4) The model works fine, suggesting that all the variables have strong, non-linear, influences on the probability of detecting the species. My main interest is in the effect dayno has on presence, given the inclusion of the other explanatory variables. Thus, I would like to extract the values of the partial prediction of dayno and its associated 2 standard errors above and below the main effect, as shown by the code plot(global_model). I have tried to extract the values by using fitted.values(global_model,dayno), but when plotted, the figure doesn't look like the partial prediction for dayno. Instead, it gives me a very scattered figure (shotgun effect). Has anyone tried to extract the partial predictions? If so, please could you advise me how to do this? Staffan P.S.. For comparison, please have a look at Simon Woods paper in R News, 1(2):20-25, June 2001, especially the figures showing the partial predictions of Mackerel egg densities. Using those graphs as an example, I would like to extract the partial predictions for e.g. s(b.depth), given the inclusion of the other variables. David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] How to extract partial predictions, package mgcv
Hi David, Yes, the strange code lines was clearly a mistake from my side. Sorry. What dataframe references did you add in your code? /Staffan David Winsemius wrote: But removing the extraneous second to last line that says just: lines ... would leave the console looking less puzzling. -- David On Sep 16, 2009, at 1:02 PM, David Winsemius wrote: It is nice to see worked examples, but there were some errors that relate to not including dataframe references. I've paste in code that runs without error on my machine: On Sep 16, 2009, at 12:39 PM, Staffan Roos wrote: Thanks Simon, I wasn't aware that predict.gam could be used in this situation. I followed you advice and managed to extract the data I needed. For people interested, I add the code with comments I used here: # # Full code for extracting partial predictions # Start the package mgcv library(mgcv) # Simulate some data (this is from Simon Wood's example in the # package mgcv manual) n-200 sig - 2 dat - gamSim(1,n=n,scale=sig) # Check the data dat ## It looks alright :-) # Run the GAM b-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat) # Plot the partial predictions (You may need to rescale your window to # see the non-linearity par(mfrow=c(1,3)) plot(b, scale=0) ### Clearly, the variables x0 and x2 are non-linear! # I would like to extract the solid line in the graphs and the # associated standard errors from the plots. Thus, I use predict # and change to a data.frame: c-data.frame(predict(b,type=terms,se.fit=TRUE)$fit) names(c)-c(pp_s.x0.,pp_s.I.x1.2..,pp_s.x2.) d-data.frame(predict(b,type=terms,se.fit=TRUE)$se.fit) names(d)-c(se_s.x0.,se_s.I.x1.2..,se_s.x2.) # Merge the three datasets: f=cbind(dat,c,d) #Restrict the number of variables to the ones I am interested in: f-f[,c(x0,pp_s.x0., se_s.x0.)] names(f) ### This data, i.e. f, can now be exported or used for further ### calculations within R #plot the data par(mfrow=c(1,1)) plot(pp_s.x0.~x0,type=p,pch=2,cex=0.50, data=f) sequence=order(f$x0) lines(f$x0[sequence],f$pp_s.x0.[sequence]) rug(jitter(f$x0)) ## Staffan, Take a look at ?predict.gam. You need to use type=terms with se=TRUE to get the quantities that you need. Simon ps. with binary data, method=REML or method=ML for the gam fit are often somewhat more reliable than the default. On Monday 14 September 2009 19:30, Staffan Roos wrote: Dear package mgcv users, I am using package mgcv to describe presence of a migratory bird species as a function of several variables, including year, day number (i.e. day-of-the-year), duration of survey, latitude and longitude. Thus, the global model is: global_model-gam(present ~ as.factor(year) + s(dayno, k=5) + s(duration, k=5) + s(x, k=5) + s(y, k=5), family = binomial(link = logit), data = presence, gamma =1.4) The model works fine, suggesting that all the variables have strong, non-linear, influences on the probability of detecting the species. My main interest is in the effect dayno has on presence, given the inclusion of the other explanatory variables. Thus, I would like to extract the values of the partial prediction of dayno and its associated 2 standard errors above and below the main effect, as shown by the code plot(global_model). I have tried to extract the values by using fitted.values(global_model,dayno), but when plotted, the figure doesn't look like the partial prediction for dayno. Instead, it gives me a very scattered figure (shotgun effect). Has anyone tried to extract the partial predictions? If so, please could you advise me how to do this? Staffan P.S.. For comparison, please have a look at Simon Woods paper in R News, 1(2):20-25, June 2001, especially the figures showing the partial predictions of Mackerel egg densities. Using those graphs as an example, I would like to extract the partial predictions for e.g. s(b.depth), given the inclusion of the other variables. David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- View this message in context: http://www.nabble.com/How-to-extract-partial-predictions%2C-package-mgcv-tp25441149p25477056.html Sent from the R help mailing list
Re: [R] How to extract partial predictions, package mgcv
Your code minus the extraneous lines was: #plot the data plot(pp_s.x0.~x0,type=p,pch=2,cex=0.50) sequence=order(x0) lines(x0[sequence],pp_s.x0.[sequence]) rug(jitter(x0)) My edition was: plot(pp_s.x0.~x0,type=p,pch=2,cex=0.50, data=f) sequence=order(f$x0) lines(f$x0[sequence],f$pp_s.x0.[sequence]) rug(jitter(f$x0)) So almost every one of the last set of plotting needed one or more additions of f$ or data=f in order to run without error. Perhaps you attached f earlier and didn't include that code? I personally have given up using attach() for just that reason. -- David. On Sep 16, 2009, at 1:31 PM, Staffan Roos wrote: Hi David, Yes, the strange code lines was clearly a mistake from my side. Sorry. What dataframe references did you add in your code? /Staffan David Winsemius wrote: But removing the extraneous second to last line that says just: lines ... would leave the console looking less puzzling. -- David On Sep 16, 2009, at 1:02 PM, David Winsemius wrote: It is nice to see worked examples, but there were some errors that relate to not including dataframe references. I've paste in code that runs without error on my machine: On Sep 16, 2009, at 12:39 PM, Staffan Roos wrote: Thanks Simon, I wasn't aware that predict.gam could be used in this situation. I followed you advice and managed to extract the data I needed. For people interested, I add the code with comments I used here: # # Full code for extracting partial predictions # Start the package mgcv library(mgcv) # Simulate some data (this is from Simon Wood's example in the # package mgcv manual) n-200 sig - 2 dat - gamSim(1,n=n,scale=sig) # Check the data dat ## It looks alright :-) # Run the GAM b-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat) # Plot the partial predictions (You may need to rescale your window to # see the non-linearity par(mfrow=c(1,3)) plot(b, scale=0) ### Clearly, the variables x0 and x2 are non-linear! # I would like to extract the solid line in the graphs and the # associated standard errors from the plots. Thus, I use predict # and change to a data.frame: c-data.frame(predict(b,type=terms,se.fit=TRUE)$fit) names(c)-c(pp_s.x0.,pp_s.I.x1.2..,pp_s.x2.) d-data.frame(predict(b,type=terms,se.fit=TRUE)$se.fit) names(d)-c(se_s.x0.,se_s.I.x1.2..,se_s.x2.) # Merge the three datasets: f=cbind(dat,c,d) #Restrict the number of variables to the ones I am interested in: f-f[,c(x0,pp_s.x0., se_s.x0.)] names(f) ### This data, i.e. f, can now be exported or used for further ### calculations within R #plot the data par(mfrow=c(1,1)) plot(pp_s.x0.~x0,type=p,pch=2,cex=0.50, data=f) sequence=order(f$x0) lines(f$x0[sequence],f$pp_s.x0.[sequence]) rug(jitter(f$x0)) ## Staffan, Take a look at ?predict.gam. You need to use type=terms with se=TRUE to get the quantities that you need. Simon ps. with binary data, method=REML or method=ML for the gam fit are often somewhat more reliable than the default. On Monday 14 September 2009 19:30, Staffan Roos wrote: Dear package mgcv users, I am using package mgcv to describe presence of a migratory bird species as a function of several variables, including year, day number (i.e. day-of-the-year), duration of survey, latitude and longitude. Thus, the global model is: global_model-gam(present ~ as.factor(year) + s(dayno, k=5) + s(duration, k=5) + s(x, k=5) + s(y, k=5), family = binomial(link = logit), data = presence, gamma =1.4) The model works fine, suggesting that all the variables have strong, non-linear, influences on the probability of detecting the species. My main interest is in the effect dayno has on presence, given the inclusion of the other explanatory variables. Thus, I would like to extract the values of the partial prediction of dayno and its associated 2 standard errors above and below the main effect, as shown by the code plot(global_model). I have tried to extract the values by using fitted.values(global_model,dayno), but when plotted, the figure doesn't look like the partial prediction for dayno. Instead, it gives me a very scattered figure (shotgun effect). Has anyone tried to extract the partial predictions? If so, please could you advise me how to do this? Staffan P.S.. For comparison, please have a look at Simon Woods paper in R News, 1(2):20-25, June 2001, especially the figures showing the partial predictions of Mackerel egg densities. Using those graphs as an example, I would like to extract the partial predictions for e.g. s(b.depth), given the inclusion of the other variables. David Winsemius, MD Heritage Laboratories West Hartford, CT __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. David
Re: [R] How to extract partial predictions, package mgcv
Yes, you're correct David, I used attach(f) earlier. Thanks for clarifying! I should change my own code accordingly. /Staffan David Winsemius wrote: Your code minus the extraneous lines was: #plot the data plot(pp_s.x0.~x0,type=p,pch=2,cex=0.50) sequence=order(x0) lines(x0[sequence],pp_s.x0.[sequence]) rug(jitter(x0)) My edition was: plot(pp_s.x0.~x0,type=p,pch=2,cex=0.50, data=f) sequence=order(f$x0) lines(f$x0[sequence],f$pp_s.x0.[sequence]) rug(jitter(f$x0)) So almost every one of the last set of plotting needed one or more additions of f$ or data=f in order to run without error. Perhaps you attached f earlier and didn't include that code? I personally have given up using attach() for just that reason. -- David. On Sep 16, 2009, at 1:31 PM, Staffan Roos wrote: Hi David, Yes, the strange code lines was clearly a mistake from my side. Sorry. What dataframe references did you add in your code? /Staffan David Winsemius wrote: But removing the extraneous second to last line that says just: lines ... would leave the console looking less puzzling. -- David On Sep 16, 2009, at 1:02 PM, David Winsemius wrote: It is nice to see worked examples, but there were some errors that relate to not including dataframe references. I've paste in code that runs without error on my machine: On Sep 16, 2009, at 12:39 PM, Staffan Roos wrote: Thanks Simon, I wasn't aware that predict.gam could be used in this situation. I followed you advice and managed to extract the data I needed. For people interested, I add the code with comments I used here: # # Full code for extracting partial predictions # Start the package mgcv library(mgcv) # Simulate some data (this is from Simon Wood's example in the # package mgcv manual) n-200 sig - 2 dat - gamSim(1,n=n,scale=sig) # Check the data dat ## It looks alright :-) # Run the GAM b-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat) # Plot the partial predictions (You may need to rescale your window to # see the non-linearity par(mfrow=c(1,3)) plot(b, scale=0) ### Clearly, the variables x0 and x2 are non-linear! # I would like to extract the solid line in the graphs and the # associated standard errors from the plots. Thus, I use predict # and change to a data.frame: c-data.frame(predict(b,type=terms,se.fit=TRUE)$fit) names(c)-c(pp_s.x0.,pp_s.I.x1.2..,pp_s.x2.) d-data.frame(predict(b,type=terms,se.fit=TRUE)$se.fit) names(d)-c(se_s.x0.,se_s.I.x1.2..,se_s.x2.) # Merge the three datasets: f=cbind(dat,c,d) #Restrict the number of variables to the ones I am interested in: f-f[,c(x0,pp_s.x0., se_s.x0.)] names(f) ### This data, i.e. f, can now be exported or used for further ### calculations within R #plot the data par(mfrow=c(1,1)) plot(pp_s.x0.~x0,type=p,pch=2,cex=0.50, data=f) sequence=order(f$x0) lines(f$x0[sequence],f$pp_s.x0.[sequence]) rug(jitter(f$x0)) ## Staffan, Take a look at ?predict.gam. You need to use type=terms with se=TRUE to get the quantities that you need. Simon ps. with binary data, method=REML or method=ML for the gam fit are often somewhat more reliable than the default. On Monday 14 September 2009 19:30, Staffan Roos wrote: Dear package mgcv users, I am using package mgcv to describe presence of a migratory bird species as a function of several variables, including year, day number (i.e. day-of-the-year), duration of survey, latitude and longitude. Thus, the global model is: global_model-gam(present ~ as.factor(year) + s(dayno, k=5) + s(duration, k=5) + s(x, k=5) + s(y, k=5), family = binomial(link = logit), data = presence, gamma =1.4) The model works fine, suggesting that all the variables have strong, non-linear, influences on the probability of detecting the species. My main interest is in the effect dayno has on presence, given the inclusion of the other explanatory variables. Thus, I would like to extract the values of the partial prediction of dayno and its associated 2 standard errors above and below the main effect, as shown by the code plot(global_model). I have tried to extract the values by using fitted.values(global_model,dayno), but when plotted, the figure doesn't look like the partial prediction for dayno. Instead, it gives me a very scattered figure (shotgun effect). Has anyone tried to extract the partial predictions? If so, please could you advise me how to do this? Staffan P.S.. For comparison, please have a look at Simon Woods paper in R News, 1(2):20-25, June 2001, especially the figures showing the partial predictions of Mackerel egg densities. Using those graphs as an example, I would like to extract the partial predictions for e.g. s(b.depth), given the inclusion of the other variables. David Winsemius, MD Heritage Laboratories West
Re: [R] How to extract partial predictions, package mgcv
Staffan, Take a look at ?predict.gam. You need to use type=terms with se=TRUE to get the quantities that you need. Simon ps. with binary data, method=REML or method=ML for the gam fit are often somewhat more reliable than the default. On Monday 14 September 2009 19:30, Staffan Roos wrote: Dear package mgcv users, I am using package mgcv to describe presence of a migratory bird species as a function of several variables, including year, day number (i.e. day-of-the-year), duration of survey, latitude and longitude. Thus, the global model is: global_model-gam(present ~ as.factor(year) + s(dayno, k=5) + s(duration, k=5) + s(x, k=5) + s(y, k=5), family = binomial(link = logit), data = presence, gamma =1.4) The model works fine, suggesting that all the variables have strong, non-linear, influences on the probability of detecting the species. My main interest is in the effect dayno has on presence, given the inclusion of the other explanatory variables. Thus, I would like to extract the values of the partial prediction of dayno and its associated 2 standard errors above and below the main effect, as shown by the code plot(global_model). I have tried to extract the values by using fitted.values(global_model,dayno), but when plotted, the figure doesn't look like the partial prediction for dayno. Instead, it gives me a very scattered figure (shotgun effect). Has anyone tried to extract the partial predictions? If so, please could you advise me how to do this? Thanks, Staffan P.S.. For comparison, please have a look at Simon Woods paper in R News, 1(2):20-25, June 2001, especially the figures showing the partial predictions of Mackerel egg densities. Using those graphs as an example, I would like to extract the partial predictions for e.g. s(b.depth), given the inclusion of the other variables. -- Staffan Roos, PhD Research Ecologist BTO Scotland __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK +44 1225 386603 www.maths.bath.ac.uk/~sw283 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] How to extract partial predictions, package mgcv
Dear package mgcv users, I am using package mgcv to describe presence of a migratory bird species as a function of several variables, including year, day number (i.e. day-of-the-year), duration of survey, latitude and longitude. Thus, the global model is: global_model-gam(present ~ as.factor(year) + s(dayno, k=5) + s(duration, k=5) + s(x, k=5) + s(y, k=5), family = binomial(link = logit), data = presence, gamma =1.4) The model works fine, suggesting that all the variables have strong, non-linear, influences on the probability of detecting the species. My main interest is in the effect dayno has on presence, given the inclusion of the other explanatory variables. Thus, I would like to extract the values of the partial prediction of dayno and its associated 2 standard errors above and below the main effect, as shown by the code plot(global_model). I have tried to extract the values by using fitted.values(global_model,dayno), but when plotted, the figure doesn't look like the partial prediction for dayno. Instead, it gives me a very scattered figure (shotgun effect). Has anyone tried to extract the partial predictions? If so, please could you advise me how to do this? Thanks, Staffan P.S.. For comparison, please have a look at Simon Woods paper in R News, 1(2):20-25, June 2001, especially the figures showing the partial predictions of Mackerel egg densities. Using those graphs as an example, I would like to extract the partial predictions for e.g. s(b.depth), given the inclusion of the other variables. -- Staffan Roos, PhD Research Ecologist BTO Scotland __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
