[R] plot in generalized additive model (GAM)

David Winsemius dwinsemius at comcast.net
Wed Jul 2 19:59:56 CEST 2014


On Jul 2, 2014, at 8:37 AM, Agostino Di Ciaula wrote:

> My aim was to perform a a generalized additive model with penalized splines to analyze the relationship between daily mortality and PM10 levels, and to control for the non-linear confounding effects of weather (temperature = “Tmax”; humidity = “umidità”).
> 
> this is my complete analysis: 
> > x.gam <- gam(mortality ~ s(PM10) + (Tmax) + (umidità), data = sab, family = quasipoisson)
> > summary (x.gam)
> 
> Family: quasipoisson 
> Link function: log 
> 
> Formula:
> mortality ~ s(PM10) + (Tmax) + (umidità)
> 
> Parametric coefficients:
>             Estimate Std. Error t value Pr(>|t|)  
> (Intercept)  1.50456    0.83802   1.795   0.0784 .
> Tmax         0.02060    0.01760   1.170   0.2472  
> umidità      0.01256    0.00826   1.520   0.1345  
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> Approximate significance of smooth terms:
>           edf Ref.df     F p-value   
> s(PM10) 3.067  3.857 4.862 0.00235 **
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> R-sq.(adj) =  0.269   Deviance explained = 34.2%
> GCV score = 6.6782  Scale est. = 5.9797    n = 58
> 
> Can I assume that the y axis of the plot (which is labeled "sPM10,3.07)
> plot(x.gam, page =1)
> is the log relative risk ?

You have repeatedly ignored my suggestions that you describe the data. Without knowing what was measured or counted in the column of data named "mortality", I doubt that anyone (and certainly not I)  can offer reassurances or advice about what to assume.

-- 
David.
> 
> thanks again
> agostino
> 
> Il giorno 02/lug/2014, alle ore 16:41, David Winsemius <dwinsemius at comcast.net> ha scritto:
> 
>> 
>> On Jul 2, 2014, at 1:25 AM, agostinodiciaula at tiscali.it wrote:
>> 
>>> Sorry for my unclear request and description
>>> I performed:
>>> 
>>> x.gam <- gam(mortality ~ s(PM10) + (Tmax) + (umidity), data = data, family = quasipoisson)
>>> plot(x.gam, page =1)
>>> 
>>> the plot showed on the left axis "sPM10,3.07".
>>> How can I have "Log relative risk" ?
>> 
>> You just want to change the axis label? Then add this to the plot call:
>> 
>> ylab="Log Relative Hazard"
>> 
>> Again, you have not told us anything about the data, so we remain unable to say whether the results corresponds to that label change, but with a quasipoisson link one might expect a logged estimate of <something>.
>> 
>> -- 
>> David
>> 
>>> thanks again for your kind help
>>> agostino
>>> 
>>> Il 02.07.2014 08:42 David Winsemius ha scritto:
>>> 
>>>> On Jul 1, 2014, at 11:33 PM, agostinodiciaula at tiscali.itwrote:
>>>>> Dear David, many thanks for your reply. The aim was to evaluate the relationship between daily mortality (variable "mortality") and mean daily values of PM10 levels (variable "PM10") in a specific geographic area, considering the effect of meteorological confounders (temperature, umidity). I used the following model (mgcv package): gam(mortality ~ (PM10) + (Tmax) + (umidity), data = data, family = quasipoisson) the question is: how can I obtain a plot of log-relative risk of daily mortality vs. PM10 levels? regards
>>>> Still rather unclear what the problem is. I see no error message in your question. Are you asking how to assign a model to a name and then call the plot method?
>>>> 
>>>> mdlAgostino
>>>> Il 02.07.2014 01:50 David Winsemius ha scritto:
>>>> 
>>>>> On Jul 1, 2014, at 12:02 AM, adc wrote:
>>>>>> I performed the following GAM by the MGCV package:
>>>>> I think it's actually spelled in all lower case.
>>>>>> gam(mortality ~ (PM10) + (Tmax) + (umidity), data = data, family = quasipoisson) How can I obtain a plot of Log-relative risk of mortality vs. PM10 ? thanks
>>>>> Shouldn't we need to know more details about the experimental setup to answer that question? And what sort of comparisons you are requesting? And about what parts of ?mgcv::plot.gam you need further explanations to answer the question?
>>>>>> agostino
>>>>> snipped
>>>>>> --
>>>>> snipped
>>>>>> Sent from the R help mailing list archive at Nabble.com. [[alternative HTML version deleted]]
>>>>> snipped
>>>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.

snipped



David Winsemius
Alameda, CA, USA



More information about the R-help mailing list