How to decompose trend and seasonality using STL in R?

Question

I have weekly data for 3 years. Now my objective is to remove the trend and seasonality effects from the series using STL function.

I can decompose time series components using decompose function in stats package. But I am getting NA values for first and last 52 values of trend and random effects.

In my sample dataset there is perfect seasonality and mean and varience are changing over time. So, I wanted to build multiplicative model. Here I have used stl function in stats package to decompose trend and seasonality. I know that stl function can handle additive model. But we can build multiplicative model also by using log transformation. Here I tried both of the models. But I am not getting as results as expected. I am sure that i am missing something in this code.

series<-ts(series,frequency=365.25/7,start(2013,9))

series<-structure(c(62L, 72L, 48L, 50L, 302L, 396L, 66L, 33L, 77L, 91L, 
38L, 38L, 43L, 45L, 134L, 754L, 1011L, 901L, 483L, 237L, 99L, 
59L, 92L, 65L, 120L, 214L, 329L, 387L, 276L, 307L, 395L, 372L, 
332L, 258L, 291L, 359L, 211L, 308L, 250L, 1374L, 1131L, 845L, 
588L, 770L, 499L, 532L, 491L, 359L, 318L, 219L, 153L, 138L, 156L, 
133L, 92L, 77L, 214L, 273L, 86L, 75L, 51L, 163L, 72L, 191L, 62L, 
49L, 79L, 573L, 569L, 444L, 410L, 404L, 345L, 141L, 146L, 179L, 
127L, 143L, 382L, 548L, 283L, 315L, 392L, 394L, 313L, 373L, 603L, 
429L, 384L, 419L, 449L, 1774L, 2025L, 1532L, 1252L, 857L, 790L, 
658L, 389L, 398L, 398L, 302L, 237L, 249L, 182L, 167L, 109L, 179L, 
377L, 288L, 146L, 126L, 449L, 138L, 580L, 130L, 94L, 150L, 173L, 
1246L, 1227L, 991L, 707L, 489L, 592L, 326L, 209L, 259L, 286L, 
243L, 344L, 335L, 368L, 397L, 349L, 313L, 1345L, 301L, 1111L, 
366L, 274L, 302L, 248L, 2518L, 2186L, 2094L, 2151L, 1847L, 1384L, 
666L, 455L, 415L, 302L, 277L, 172L, 186L), .Tsp = c(1, 3.97056810403833, 
52.1785714285714), class = "ts")


#Model 1

model1<-stl(series,"periodic",robust="TRUE")

op<-as.data.frame(model1$time.series)

head(op,25)

matplot(op,type="l")

#Model 2

model2<-stl(log(series),"periodic",robust="TRUE")

op<-exp(as.data.frame(model2$time.series))    

matplot(op,type="l")

How can I improve the model performance?

Please suggest me if there are any better ways to solve with this problem.

Thanks in advance.