[R] plot in generalized additive model (GAM)

[David Winsemius]

On Jul 2, 2014, at 1:28 PM, Agostino Di Ciaula wrote:

> "mortality" is the daily number of deaths.
> many thanks again for your patience and for your help.
> agostino

A risk estimate requires specification of a risk set, which I do not see in the data description. Usually one would add an offset of the log(risk_count). It is possible to estimate relative risks without defining an explicit risk set under certain conditions, but it is not yet clear whether those conditions would be met for whatever data collection method or design was used.

-- 
David.


> 
> 
> Il giorno 02/lug/2014, alle ore 19:59, David Winsemius <dwinsemius at comcast.net> ha scritto:
> 
>> 
>> 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
> 

David Winsemius
Alameda, CA, USA