title: Intercept calculation in Theil-Sen estimator

text: I've just written some C code for Theil-Sen, after some Googling (I don't have any definitive documentation on it). My understanding of the intercept calculation is that I first calculate the median slope, and then construct a line through every data point with this slope, find the intercept of every line, and then take the median intercept. The only way I can find to test the code is to compare the results against the Kendall-Theil Robust Line program, from the USGS. On a dataset of 237 points (healthcare data, with a Pearson correlation of ~0.55), we agree exactly on the median slope, but disagree on the intercept (by 1.4%). According to my figures, the KTRL intercept isn't the median intercept, but is instead 46% of the way through the range. After some digging around in the KTRL code, it appears that they calculate the intercept by creating a single 'median line', rather than the median of all intercepts. Their intercept is `medianY - medianX * median slope`. Any feedback on which is the "right" way to do this, if there is one, or how this is handled in R/etc? Thanks.


title: Mixing User Data For Cross-Validation

text: I have data from 12 users, and want to perform cross-validation. Is it necessary that I create my training and test data from different users, or can all of the data be randomized, and then split into Train and Test? **EDIT**   This dataset is large, as we are observing user's typing patterns, for multiple typing sessions. We're measuring 100s of features, once per answer per user. Hopefully this helps give a better idea of our data.


title: Difference between a census and a sample

text: I know when you do a census, you cover the entire population but when you take a sample, you select a group from the population. There can be only one census but there can be many different samples from the same population. Hence, the latter results in sampling error (because it is just one aspect of the population). This may be a bit simplistic but I think it captures the essential differences between a census and a sample. A consultant that my company has engaged is blurring this distinction between a census and a sample. He has essentially collected a sample, weighted it and then interpreted the findings as if they relate to the entire population. He is mostly using descriptive statistics. He has even calculated the sampling error. My question is: Is it ever valid to treat a sample as a census? If so, what are the preconditions?


title: How to read the the goodness of fit on nls of R?

text: I am trying to interpreting the output of nls(). I have read this post but I still don't understand how to choose the best fit. From my fits I have two outputs:               > summary(m)            Formula: y ~ I(a * x^b)            Parameters:       Estimate Std. Error t value Pr(>|t|)           a 479.92903   62.96371   7.622 0.000618 ***       b   0.27553    0.04534   6.077 0.001744 **        ---       Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1             Residual standard error: 120.1 on 5 degrees of freedom            Number of iterations to convergence: 10        Achieved convergence tolerance: 6.315e-06       and               > summary(m1)            Formula: y ~ I(a * log(x))            Parameters:       Estimate Std. Error t value Pr(>|t|)           a   384.49      50.29   7.645 0.000261 ***       ---       Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1             Residual standard error: 297.4 on 6 degrees of freedom            Number of iterations to convergence: 1        Achieved convergence tolerance: 1.280e-11      The first one has two parameters and smaller residual error. The second only one parameter but worst residual error. Which is the best fit?


title: What will be/how do I estimate US household income in 2020?

text: I would like to know what US household income will be in 2020 by age group (eg. 15-34,35-44,45-54,65+). The US Census provides a historical Excel sheet (http://www.census.gov/hhes/www/income/data/historical/household/2010/H10AR_2010.xls) that breaks down this data over the last few decades but of course stops at the year 2010. I've searched the web wide and far and can't seem to find projections of any kind. Is there a way I can calculate the 2020 figures myself (obviously resulting in a very rough estimate) or is there a projection somewhere out there? Thanks!


title: Why do I need bag composition to calculate OOB error of combined random forest model?

text: Could someone please explain me the answer of the question already given here: Combining randomForests in R, why are the err.rate, mse and rsq components NULL I do not understand the phrase "clip it to out-of-bag" from model$inbag.


title: Simple algorithm for online outlier detection of a generic time series

text: I am working with a large amount of time series. These time series are basically network measurements coming every 10 minutes, and some of them are periodic (i.e. the bandwidth), while some other aren't (i.e. the amount of routing traffic). I would like a simple algorithm for doing an online "outlier detection". Basically, I want to keep in memory (or on disk) the whole historical data for each time series, and I want to detect any outlier in a live scenario (each time a new sample is captured). What is the best way to achieve these results? I'm currently using a moving average in order to remove some noise, but then what next? Simple things like standard deviation, mad, ... against the whole data set doesn't work well (I can't assume the time series are stationary), and I would like something more "accurate", ideally a black box like: double outlier_detection(double* vector, double value); where vector is the array of double containing the historical data, and the return value is the anomaly score for the new sample "value" .


title: Getting rid of sparks in sample data

text: > **Possible Duplicate:**   >  Simple algorithm for online outlier detection of a generic time series   >  Getting rid of spikes in sample data How could I get rid of sparky (aka spikey) data in a discrete data set, but in a "smoother out" manner? Take for instance ![enter image description here](http://i.stack.imgur.com/vb44Q.png) There are two sparks, at 20000, but the next one at 600 is also considered a spark. I've managed to get the very high ones to zero, by               a = 2     b = 5     beta_dist = RealDistribution('beta', [a, b])     f(x) = x / 19968     normalized_insertions = [f(i) for i in insertions]          insertions_pairs = [(i, beta_dist.distribution_function(i)) for i in normalized_insertions]     plot_b = beta_dist.plot()          show(list_plot(insertions_pairs)+plot_b)      No idea how to go about the lower ones. The maximum should be reached at 100, perhaps the parameters for the beta distribution need a little more twiddling? Currently, it looks like this: ![enter image description here](http://i.stack.imgur.com/7hKHT.png) If possible, use sage as a reference for your explanations.


title: Finding anomalies using moving average in a time series

text: I want to find anomalies in a time series. Is it possible to find anomalies using moving average?


title: Spotting anomalies in time series

text: > **Possible Duplicate:**   >  Simple algorithm for online outlier detection of a generic time series I've got data on quantities sold and average price, by date, for a number of commodities and a number of outlets. Because of the provenance of the data, some of the observations are wrong. I'd like to be able to flag the potential anomalies so that they can be checked. I've done some context-independent analysis (looking for values outside specified bounds), but I'd like to try some context-dependent analysis - for example, price for a particular commodity out of line by looking across outlets for a single date, or quantity out of line by looking across dates for a single outlet. For the latter, I'm thinking that some sort of moving window would be the thing to use - I'm pretty sure that for some commodities there are both long term trends and seasonal variation going on, so what I'm looking for is spikes. My current thoughts are something along the lines of: for each observation, flag it if it's more than some distance from the mean or median of the comparator observations, where the comparator observations would be either those for the same commodity and outlet in a time window, or those for the same date and commodity for all (other) outlets. Is this anything like any standard technique? If so, what's its name, and what drawbacks does it have? (Once I've got a name I'll be able to search for more descriptions, and hopefully an R implementation). Thanks for any pointers.


title: Application to a real life problem: identifying outliers in univariate time series data

text: I have a rather simple problem that I'm having trouble deciding an answer upon. As a student studying Statistics I'm very familiar with terminology and theory but I suppose I'm stumped on my first real application of classroom teachings. I have data on the size of database uploads that occur everyday. The database is ever increasing with no upper bound. Right now I'm working with about 2 months of information and we'll assume I want to keep collecting records indefinitely. eg: 03/06/13 - 239238 04/06/13 - 240594 05/06/13 - 256920 etc... I want to generate an alert if either 1) the current day's file size is smaller than the previous day's file size 2) the current day's file size is unusually larger than expected compared to the previous day's file size I want to stress that I merely only want an alert, or to know of when either case occurs, not to delete the record of this instance or do anything to it, as many of the searches I've tried to find of an existing question similar to mine has to do with how to handle the outlier itself. So the problem for me is, what measure should be used to more rigorously determine if today's file size is an outlier as opposed to some arbitrary rule? I ideally want to use the given historical data rather than just saying 'if today's size is 10% bigger than yesterdays, alert me'. Is the boxplot good enough? Is calculating and using the mean and standard deviation 'better' (browsing around would leave me to believe that the latter is outdated and the former is better)? It is entirely possible that the database would need to be cleaned leading to a large drop, or that the database would have a huge import leading to a large increase, both of which are acceptable but would throw off the mean and standard deviation. Any help is appreciated. Thank you.


title: Getting rid of spikes in sample data

text: > **Possible Duplicate:**   >  Simple algorithm for online outlier detection of a generic time series How could I get rid of sparky data in a descrete data set, but in a "smoother out" manner? Take for instance http://i.stack.imgur.com/vb44Q.png There are two sparks, at 20000, but the next one at 600 is also considered a spark. I've managed to get the very high ones to zero, by               a = 2     b = 5     beta_dist = RealDistribution('beta', [a, b])     f(x) = x / 19968     normalized_insertions = [f(i) for i in insertions]          insertions_pairs = [(i, beta_dist.distribution_function(i)) for i in normalized_insertions]     plot_b = beta_dist.plot()          show(list_plot(insertions_pairs)+plot_b)      No idea how to go about the lower ones. The maximul should be reached at 100, perhaps the parameters for the beta distribution need a little more twiddling? Currently, it looks like this: http://i.stack.imgur.com/7hKHT.png If possible, use sage as a reference for your explanations.


title: Specifying a model in MCMCglmm

text: I am trying to fit a GLMM: logit$[P(Y_{ij}=1)] = \beta_0 + \beta_1 X_{ij} + \beta_2R_i + \gamma_i$ where $\gamma_i\sim N(0,\sigma^2_\gamma)$, $j=0,1$ and $i=0,1,...,94$. I first went with the lme4 package in R:               m <-  glmer(PHENO==2 ~ GENO + RACE + (1|GROUP), family="binomial",data=iih_data,nAGQ=1 )      but this uses the Wald z-test by default for hypothesis test for $\beta_1=0$ which is probably not applicable in this case. It was suggested here http://glmm.wikidot.com/faq to use either MCMC or parametric bootstrap to obtain more valid p-values. I rolled my own bootstrap procedure, but it was rather slow. So I wanted to see if MCMC method would be faster. This led me to the R package MCMCglmm. However, this package is above my understanding of statistics and I was unable to follow the examples on how to set up this analysis using MCMCglmm. Can you help me set up the above analysis using MCMCglmm to obtain a valid p-value for $\beta_1$?


title: Estimating the autocorrelation function with some noisy observations

text: I am wondering how to estimate the actual correlation function when I have some noisy samples of some space. Lets say, I have a space and the locations in the space are variables following a multivariate gaussian distribution. I generate some samples with the given autocorrelation function and add some noise independently to all of them. Now from the samples I find an empirical correlation function and then fit a model autocorrelation function to it. This model might has some nugget effects or might not represent the true autocorrelation function. How to deal with this issue? For eg, I generated some data with exponential correlation function with range 10. If I add no noise and try to fit a model correlation function(exponential) based upon some samples of the field. I get almost the actual function. I generated data with exponential correlation function with range 10. I then got some samples from this function. I got the empirical correlation function and then fit a model (exponential) to it to get almost the same params i.e range of 9.5 However, when I add some white gaussian noise to the correlation/covariance function(I assume std to be 1) and then I tried to fit a model correlation function to the empirical one, I got range of 5 and also nugget of 0.8. How to estimate the underlying function in this case?


title: Posterior probability - change in Beta hyperparameters

text: Can you explain, how does $\text{B}(\alpha, \beta)$ transfrom to $\text{B}(s+\alpha, f+\beta)$ in the following equation? $$ \begin{align*} p(\left. q=x \right| s,f) &= {{{s+f \choose s} x^{s+\alpha-1}(1-x)^{f+\beta-1} / \operatorname{B}(\alpha,\beta)} \over \int_{y=0}^1 \left({s+f \choose s} y^{s+\alpha-1}(1-y)^{f+\beta-1} / \operatorname{B}(\alpha,\beta)\right) dy} \\\ ~\\\ ~\\\ ~\\\ ~\\\ & = {x^{s+\alpha-1}(1-x)^{f+\beta-1} \over \operatorname{B}(s+\alpha,f+\beta)} \end{align*} $$ (This equation is taken from Wikipedia's article on conjugate priors.) However, I started an investigation from parameter estimation for text analysis by G. Heinrich (in that paper, equations 23 and 24). I see, that in many papers, this is an "of course" statement, but I would be grateful for a derivation of the formula, or hints on how to do it on my own.


title: How to calculate SE for the ratio of two proportions?

text: I have evaluated the prevalence (proportion) of a disease in men and women, in different studies available in the literature. So I have sample sizes for men and women (in different studies) and also know the number of the affected people in each gender. So I can calculate standard errors for each prevalence. However, the more interesting finding is the ratio of the prevalence (proportion) of disease in males divided by its prevalence (proportion) in females (male:female ratio). I have calculated these too. However I need the SE for those "ratios of proportions" in order to draw a funnel plot. I couldn't find a proper way to calculate SE for those ratios. Do you think it is possible? If not, how am I supposed to draw funnel plots for my male:female ratios? Many thanks.


title: Pre-post analysis with multiple post measures

text: I am having trouble deciding what type of analysis to use for the following data: I have two groups (one control, one treatment) and took a pre-measure of my dependent variable, and then 3 post-measures. I am interested in understanding if the treatment affected the dependent variable over the course of the study. I originally did a rmANOVA, and found no significant interaction (treatmentxtime), but found a main effect of the treatment. I think this means that the treatment group was different than the control over the entire study, including before treatment began. My question is, is there a better analysis that takes into account if the groups were in fact different before treatment? With a t-test, they were not different, but somehow with all of the data throughout the experiment, they are different. Thanks!


title: Kolmogorov distribution

text: Is there a package in R for the Kolmogorov distribution which allows me to plot density, distribution, calculates quantiles, etc.? The Kolmogorov distribution arises from $K=\sup|B|$, where $B$ is a Brownian bridge. Its values are usually tabulated, so I thought it would have its own function in R, like the normal distribution. It seems `ks.test()` doesn't use a R function for the quantiles, but actually calls a C++ function.


title: Statistical methods to compare spatial maps

text: I would like to compare 4 radon flux maps (Bq m-2 s-1) generated using different methods. All of the maps use same projection and have same resolution. I am having trouble finding statistical method(s) which i could use for comparing the maps and analyse regions where they agree well and regions where the difference is large. I would be thankful for any ideas/suggestions that I can have a look at for further study. My final goal is to merge all these maps together and create a new map and run some simulations for the final map. Can i just take mean of the respective cell values to make the new map. Is it a correct way? Thanks a lot for your time and help. And off course, a reproducible example:               library(raster)               e = extent(1, 1+(5*1), 1, 1 + (5*1))      r1 = raster(e)      res(r1)=c(1,1)      set.seed(1)     r1[] = runif(5*5)      plot(r1, main='map1')           e = extent(1, 1+(5*1), 1, 1 + (5*1))      r2 = raster(e)      res(r2)=c(1,1)      set.seed(2)     r2[] = runif(5*5)      plot(r2, main='map2')          e = extent(1, 1+(5*1), 1, 1 + (5*1))      r3 = raster(e)      res(r3)=c(1,1)      set.seed(3)     r3[] = runif(5*5)      plot(r3, main='map3')               e = extent(1, 1+(5*1), 1, 1 + (5*1))      r4 = raster(e)      res(r4)=c(1,1)      set.seed(4)     r4[] = runif(5*5)      plot(r4, main='map4')      ![enter image description here](http://i.stack.imgur.com/hUelA.jpg)


title: How to compute K and n? [Item-based Collaborative Filtering]

text: I'm currently studying an item-based collaborative filtering algorithm described in Ul Haq, Raza - Hybrid Recommender System Towards User Satisfaction. I've formulated the algorithm below based on it. I have no problem on steps 1 to 3 but ...   * in step 4 it says there: **Set the threshold value n.** How can I determine the value of n? Is there a formula for getting the value of it? I already checked it out but there's nothing there.    * in Step 5: **Select the K most similar products in M**. The same question, how can I compute the value of K? >   1. Retrieve all the item rated by an active user and put it to Q. >   2. Isolate the users who have rated both the target item (i) and the items > rated by the active user in Q, get the item and put it in R. (co-rated > items) >   3. Calculate the item similarities using the Pearson Correlation > Coefficient with all the items (j) in R. >   4. Set the threshold value n, If the similarity of i and j is greater or > equal to n, (sim(i,j) >= n), Then include it in M. >   5. Select the K most similar products in M. >   6. Take the weighted average of the users rating on these similar items K. >


title: Sample variance robustness

text: I am trying to understand the robustness of the sample variance. I want to calculate its influence function and in order to do so a previous step is to obtain the functional for the contaminated distribution: $$T_{\sigma^2}(F_\epsilon) =T_{\sigma^2}((1-\epsilon)F+ \epsilon \delta_x)$$ In particular, there are two things I do not fully understand regarding the derivation of $T_{\sigma^2}(F_\epsilon)$:   1. The references I am using state that the sample variance functional can be expressed as: $$T_{\sigma^2}(F) =\int{(y-T_{\mu}(F))^2dF(y)}= \frac{1}{2} \int\int (y-z)^2dF(y)dF(z)$$ Could anyone explain me how the second equality is obtained?   2. Using the previous expression, the functional $T_{\sigma^2}(F_\epsilon)$ can be obtained as: $$T_{\sigma^2}(F_\epsilon) =T_{\sigma^2}((1-\epsilon)F+ \epsilon \delta_x) = (1-\epsilon)^2 T_{\sigma^2}(F) + \epsilon(1-\epsilon) \int (y-x)^2dF(y)$$ In this case, I do not fully understand where the term $\epsilon(1-\epsilon) \int (y-x)^2dF(y)$ comes from. Once these equalities have been obtained, I have no further problems to derive the influence function for the sample variance, which final result is $IF(x, T, F) = (x-\mu)^{2} - \sigma^2$. **EDIT:** two small typos corrected


title: Testing normality assumptions for linear mixed models and mixed (repeated) GLM ANOVA in SPSS

text: I have a mixed design that includes both repeated (condition) and between (sex and genotype) subjects factors. I would like to assess whether my data meets the normality assumptions for 1) General linear models (repeated) and 2) linear mixed models using SPSS. What is the best method for doing so? Assuming that I find violations to the normality assumptions how do I select the optimal method(s) for transformation?


title: Expected sum of money

text: There is a little problem in the working on this one I have encountered recently. Here it goes: Starting off with \$100, when you toss a coin, if there is a tail, you will get \$10 and lose $10 if it is a head. What is the probability that after 1000 tosses, the amount of money you have is between 90 and 110? Here is the way I do it. Let $X$ be the difference in amount when you toss a coin. Hence, X could be -10 and 10 with the probability of 1/2 each. The expected value is 0 and the standard deviation is 10. Hence, I go on formulating the sum $S = X_1 + X_2 + \ldots + X_{1000}$ and figure out the $\mu(S) = 0$ and $var(S) = 100000$ and I do not know how to go from here? Any hints will be greatly appreciated.


title: R one vs all classification example code

text: I'm trying to build a one vs all multi-class classifier in R and my results are pretty poor. I'm curious if anyone can suggest a piece of example code somewhere on the internet. I'm trying to classify some text. There are 5 classes that the text strings could fall into. When I search for 'one vs all classifier' I get pages of papers on how great they are but can't find a single implementation.


title: Retrospective case study and inferential statistics

text: I am carrying out some research and am a novice to statistics, having done a 2 week module primarily focusing on parametric data and experimental design. However, for my research I am doing a retrospective case study, on the hearing and speech outcomes of 3 groups of patients with 3 different ear conditions following a surgical intervention to improve their hearing. Ear condition being the independent variable, speech and hearing outcomes, the dependant variable. The sample size is based on available data not a power calculation and each group has the following numbers, each group represents a different ear condition: 1) n =10, 2) n= 5 and 3) n=2. I have been advised to do inferential stats, so have done the appropriate tests, but I wonder if inferential stats are geared toward an experimental design? The outcomes are also affected by many other demographic factors between patients such as the age they had the surgery, underlying cognitive ability and so forth- these were not controlled for owing to the type of study. I also wonder as no power was calculated before (normally it is done following a pilot) but would you do a post hoc power calculation in this type of study? I have also done descriptive stats, but I am quite uncertain about the use of inferential stats (Mann Whitney and t-tests) were appropriate to the data...but I am unsure what inferences can be made on this type of study design? Any thoughts would be greatly appreciated.


title: EM algorithm with a component for outliers

text: i have a vector of measurements from one to three classes, which can be modeled by gaussian distributions. There are some outliers in the data. I use the EM algorithm to learn the parameters of the components. One property of the algorithm is, that the incomplete likelihood of the data increases monotonically. My code works fine, but the parameters are heavily biased due to the outliers. Now i want an extra component, a uniform distribution, to model the outliers. And now the monotonicity is no longer guaranteed. As far as i know the new component only changes the E step of the algorithm. Does anyone know what my mistake is?               library(gtools)          #create 200 data points from two components plus one outlier     y = c(rnorm(n=100,mean=-4,sd=.5),rnorm(n=100,mean=0,sd=.5),4)     stripchart(y,method='jitter',pch=4)               em = function(x,maxit){          K = 2 #number of components     n = length(x) #size of dataset          #start values     p = 1 / (max(x)-min(x)) #probability for the uniform distribution     pi = rdirichlet(n=1,alpha=c(rep(5,times=K),1))[1,] # mixing probabilities     mu = sample(x,size=K) #vector of means     sigma = .7 #standard deviation     gamma = 0 #           #maximization     m_step = function(){         pi = apply(gamma,2,sum)         mu = gamma[,1:K] * x         mu = apply(mu,2,sum)         mu = mu / pi[1:K]         pi = pi / n         sigma = t(sapply(x,`-`,mu))         sigma = gamma[,1:K] * sigma^2         sigma = sqrt(sum(sigma) / n)         pi <<- pi; mu <<- mu; sigma <<- sigma         return()     }          #expectation     e_step = function(){         gamma = sapply(1:K,function(k){ dnorm(x,mean=mu[k],sd=sigma)},USE.NAMES=FALSE)         gamma = cbind(gamma,p)         gamma = t(t(gamma)*pi)         tmp = apply(gamma,1,sum)         gamma <<- gamma / tmp         loglik = sum(log(tmp))     }          loglik = rep(NA,times=maxit)          loglik[1] = e_step()     m_step()          j = 1     gain = Inf     while(j <= maxit && gain > 1e-8){         loglik[j+1] = e_step()         gain = loglik[j+1] - loglik[j]         m_step()         j = j+1     }          if(any(diff(loglik[1:j])<0)){         print('Error: log likelihood not monotonic')         print(loglik[1:j])     }          return(list(pi=pi,mu=mu,sigma=sigma,loglik=loglik[j]))     }          print(em(y,200))      I found a solution. I introduced a new hidden binary variable with the meaning outlier/no outlier and now it works. I will try to write down the math and post it later.


title: Simplified table for computing higher order cumulants

text: I'm looking for a simplified table that shows how to compute the cumulant cum(x1, x2, ... xn) for **multiple** variables, n > 4. I know that the 2nd and 3rd order cumulants are equivalent to the compound moments of the same order. For the 4th order cumulant one also has to include the partitions such that after simplification one gets: cum(x1,x2,x3,x4) = E{x1,x2,x3,x4} - E{x1,x2}E{x3,x4} \- E{x1,x3} * E{x2,x4} \- E{x1,x4} * E{x2,x3} according to: http://sipi.usc.edu/~mendel/publications/HOSATutorial.pdf (This source also shows how you can generically get the cumulant for higher orders (p.297) but just applying this term would result in very long terms, making it infeasible to use it in real computations). I am hoping that such a table already exists somewhere, given that there are also other people who do work with higher order cumulants.


title: plotting part of data using violin plot and the other part as stripchart (stripplot)

text: I an trying to plot 6 groups, in which 3 has a good amount of counts (ranging from 20 to 400) and the rest 3 have very small counts (1 to 9). I want to plot the distribution of a continuous variable within each of these 6 groups and be very honest in the plot about the dramatically different sizes of the groups. I am thinking about generating violin plot for the 3 groups with larger counts to show both the median and the dispersion, and for the low-count groups, plot stripchart to show individual data points. Could anyone kindly help me find a way to do such a hybrid plot? Other plot suggestions are also welcome.


title: Confirmatory factor analysis about technology acceptance model in R

text: I tried to do a cfa with lavaan package in R. I'm not sure, if SO is the right stackexchange plattform, if mathematics fits better just tell me and i will move it. Here is my model: ![enter image description here](http://i.stack.imgur.com/2dg6h.jpg) I know how to define my Latent variables like this:               > MASTER.model <- ' Usability  =~ PUS1 + PUS2 + PUS3 + PUS4 + PUS5      +                   PerceivedUsefulness =~ PU1 + PU2 + PU3 + PU4 + PU5     +                   PerceivedEaseOfUse   =~ PEOU1 + PEOU2 + PEOU3 + PEOU4 + PEOU5     +                   BehavioralIntent =~ BI1 + BI2      But I don't know how to define the links between Usability, EOU, etc. Is this a covariance, written with e.g. `Usability~~PerceivedEaseOfUse` or do I write it like this: `PerceivedEaseOfUse~Usability` which means PerceivedEaseOfUse is regressed from Usability? My try was:               MASTER.model <- ' Usability  =~ PUS1 + PUS2 + PUS3 + PUS4 + PUS5      +                PerceivedUsefulness =~ PU1 + PU2 + PU3 + PU4 + PU5     +                   PerceivedEaseOfUse   =~ PEOU1 + PEOU2 + PEOU3 + PEOU4 + PEOU5     +                   BehavioralIntent =~ BI1 + BI2     +                   PerceivedUsefulness ~ PerceivedEaseOfUse     +                   BehavioralIntent ~ PerceivedEaseOfUse     +                   BehavioralIntent ~ PerceivedUsefulness     +                    Usability ~~ Usability     +                    PerceivedEaseOfUse ~~ PerceivedEaseOfUse     +                    PerceivedUsefulness ~~ PerceivedUsefulness                      +                    PerceivedEaseOfUse ~ Usability'      This gave me a results, but my RMSEA was realy poor which would mean my model doesn't fit my data. I just want to be sure it is not wrong because of me using R in a wrong way. I would like to get results similar to these in the following picture: ![enter image description here](http://i.stack.imgur.com/DbyVu.png) Thank you.


title: Variables in factor analysis with different metrics/response categories

text: If my data has different response categories (i.e. a 3-category vs. 4-category ordinal data), can I use standard factor analysis procedures? Do I need to transform my data in any way before starting analysis? Thank you!


title: Simple and paired t-tests and the control for multiple comparisons

text: Let's say that I have a group of subjects that were measured in two conditions (A and B) and I want to check three hypotheses: i) do data from A and B differ from each other?, ii) do data from A differ from zero?, iii) do data from B differ from zero?. So, my tests would be: a paired t-test for hypothesis i), and two simple t-tests for hypotheses ii) and iii). I'm wondering whether I should or I shouldn't correct the alpha level of these three tests. I'm inclined to apply the correction for multiple comparisons to avoid type I errors, but at the same time I would also like to avoid type II errors.


title: Using R to bootstrap the variance difference

text: I asked this question on Stack Exchange, but I think it might be too specialized. Hopefully someone in the mixed model group can help me out. I want to be able to bootstrap the variance differences between two data sets obtained at different times while taking out the error in a random effect. I have 2 sets of experimental data, where the data was measured at 2 time points (initial and final). I also have a set of simulation data. I want to compare the variance of the simulated date with the variance difference between the experimental data (final - initial). The idea is to get confidence intervals from the bootstrap to compare the experimental data with the simulation. I am having trouble making the statistic for the bootstrap function in the `boot` package for R. So far I have.               varcomp <- function ( formula, data, indices ) {         d <- data[indices,] #sample for boot         fit <- lmer(formula, data=d) #linear model         res.var = (attr (VarCorr(fit), "sc")^2) # variance estimation         return(res.var)         }      But this function only returns the variance of a single data set. I want to be able to input 2 sets of data and have it return the difference between the two data sets' variance. When I try something like:               varcomp <- function ( formula, data1, data2, indices ) {     d1 <- data1[indices,] #sample for boot     d2 <- data2[indices,] #sample for boot     fit1 <- lmer(formula, data=d1) #linear model     fit2 <- lmer(formula, data=d2) #linear model     a = (attr (VarCorr(fit1), "sc")^2) #output variance estimation     b = (attr (VarCorr(fit2), "sc")^2) #output variance estimation     drv = a - b #difference between the variance estimations     return(drv)     }      I would then put it into boot such as:               ip1.boot <- boot ( data = ip1, statistic=varcomp, R=100, formula=CNPC~(1|Cell.line:DNA.extract)+Cell.line)      I can't do it this way because the boot function only allows for one data set to be inputted. **Does anyone know how to create the correct statistic function for this?** An example of the data can also be downloaded here (2 csv files zipped 1.22KB.) My data looks something like the following: Initial                      Cell.line    Time DNA.extract   Gene      CNPC     1          9 initial           1 atubP1 1778.4589     2          9 initial           1 atubP1 2108.0552     3          9 initial           1 atubP1 2118.6725     4          9 initial           2 atubP1 2018.6593     5          9 initial           2 atubP1 1935.9008     6          9 initial           2 atubP1 1749.9158     7          9 initial           3 atubP1 1524.7475     8          9 initial           3 atubP1 1532.9781     9          9 initial           3 atubP1 1693.3098     10        17 initial           1 atubP1 1076.4720     11        17 initial           1 atubP1 1101.3315     12        17 initial           1 atubP1 1185.3606     13        17 initial           2 atubP1 1131.1118     14        17 initial           2 atubP1  892.7087     15        17 initial           2 atubP1 1028.5465     16        17 initial           3 atubP1  887.9972     17        17 initial           3 atubP1  732.9646     18        17 initial           3 atubP1  680.6724      Final                  Cell.line  Time DNA.extract   Gene      CNPC     1          9 final           1 atubP1 1262.2378     2          9 final           1 atubP1 1261.9858     3          9 final           1 atubP1 1390.6873     4          9 final           2 atubP1 1539.7180     5          9 final           2 atubP1 1510.5405     6          9 final           2 atubP1 1443.1767     7          9 final           3 atubP1 1456.2050     8          9 final           3 atubP1 1578.6396     9          9 final           3 atubP1 1656.1822     10        17 final           1 atubP1 1462.5179     11        17 final           1 atubP1 1580.9956     12        17 final           1 atubP1 1255.9020     13        17 final           2 atubP1  886.7579     14        17 final           2 atubP1  581.8116     15        17 final           2 atubP1  722.0526     16        17 final           3 atubP1 4168.7895     17        17 final           3 atubP1 3266.2105     18        17 final           3 atubP1 4219.5645


title: Prediction errors for Adaboost

text: I understand how the adaboost algorithm works to produce a prediction of a class, however one thing I haven't seen is how to get a measure of accuracy for that prediction. For example, if I fit a logistic regression, then I get a set of coefficient estimates _and their standard errors_. This allows me to create a prediction interval for a new unseen case, as well as a point estimate. However, for adaboost, I've only ever seen formulae and algorithms for creating the prediction, but nothing for the prediction interval. The only thing I can think of is the test set error - but this is a "global" error estimate. Another, more mathematical way to formulate this question is: Suppose we fit an adaboost algorithm to data $D_{T}\equiv\\{(x_1,y_1),\dots,(x_T,y_T)\\}$ and get an estimated prediction function $\hat{f}$. Now we use this to predict $n$ new observations $(\hat{f}(x_{T+1}),\dots,\hat{f}(x_{T+n}))$ (we don't observe $y$ for these data points). What is the mean square error matrix $MSE(\hat{y})=E((y-\hat{y})(y-\hat{y})^T)$?


title: How to choose initial values for nonlinear least squares fit

text: The question above says it all. Basically my question is for a generic fitting function (could be arbitrarily complicated) which will be nonlinear in the parameters I am trying to estimate, how does one choose the initial values? I am trying to do nonlinear least squares fit by the way. Is there any strategy or method? Has this been studied? Any references? Anything besides ad hoc guessing? Specifically, right now one of the fitting forms I am working with is a Gaussian plus linear form with five parameters I am trying to estimate, like $$y=A e^{-\left(\frac{x-B}{C}\right)^2}+Dx+E$$ where $x = \log_{10}$(abscissa data) and $y = \log_{10}$(ordinate data) meaning that in log-log space my data looks like a straight line plus a bump which I am approximating by a gaussian. I have no theory, nothing to guide me about how to initialize except perhaps graphing and eyeballing like the slope of the line and what the center/width of the bump is. But I have more than hundred of these fits to do so instead of graphing and guessing, I would prefer some approach which can be programmed and executed by my machine. The only thing I can think of is to just randomly choose initial values. MATLAB offers to choose values randomly from [0,1] uniformly distributed. So with each data set, I run the randomly initialized fit a thousand times and then pick the one with the highest $r^2$. Any other (better) ideas? Looking for something as generic as possible, if not then maybe for my fitting form above and some simple curves like power law, etc. Another not-so-unrelated question, how do you describe your "normal" space where the data was originally collected as opposed to log-log space? So after the log-transform, you are in log-space. Where are you before the log- transform? Basically, how do you fill in the blank, "In **___ __ ___ ___** space we have a perfect power law which in log-space is a straight line". Thank you.


title: Cumulative Incidence vs. Kaplan Meier to estimate probability of failure

text: To estimate the probability failure in medical sciences, it is not atypical to use 1-KM. However this does not account for competing risks, such as death by natural causes or causes unrelated with the disease, which preclude the event of interest. Thus 1-KM provides an inadequate measure, and cumulative incidence curves, such as the ones used in cmprsk (in R). My question is that a lot of medical literature still reports KM curves or even 1-KM. Does this mean that the results reported in a lot of medical literature may be inaccurate (or more precisely over-estimated)? Or are there reasons why 1-KM would be preferred? Furthermore, if there is a difference between 1-KM and the cumulative indices curve what other parts of your analysis are also effected (i.e. discrimination, calibration...)?


title: Mean or median as the measure of central tendency for data spanning three orders of magnitude

text: Here is a data set I am working with. ![enter image description here](http://i.stack.imgur.com/Smxen.png) I basically have a total of 1082 different data segments of $y$ vs. $x$ all plotted here on a log scale in bright blue in the background. Now what I need to do is to choose one "representative" data segment from all these data segments. Each data segment consists of 20 points. For comparison, I have computed and plotted the $x$-wise (meaning for each $x$, consider all 1082 $y$) arithmetic mean (in magenta), geometric mean (in yellow), and the median (in black) as well as the interquartile range. The problem is that for any given $x$, $y$ values span about three orders of magnitude and my intuition tells me that the arithmetic mean doesn't make sense here. We obviously have some outliers (both above and below) but the higher outliers are very large so they severely skew the arithmetic mean. My instinct says to choose the median. My question is, what is the appropriate measure of central tendency here? How should it be computed? I would like a better reason than just my intuition. Is the median appropriate here? Maybe an entirely different measure of central tendency than these three? I understand that my "central" data segment will most likely may NOT be one of the actual data segments which is okay. If there are any references which discuss this, they would be appreciated too. Thank you.


title: Analysis of Likert scale

text: > **Possible Duplicate:**   >  Analysis of part of Likert scale My survey consist of 6 airlines (non metric) and 10 attributes for rating these airlines. The questionnaire uses 5 point likert scale. Only one airline can be rated using one questionnaire. The questionnaire has 10 questions, each one regarding one attributes. My concerns are:   * What are the dependent and independent variable?   * Whether the attributes rating using Likert scale are metric values or non metric values.   * My final objective is to sort out the best airline based on the ratings obtained from the passengers. Based on the dependent and independent variables which analytical tool have to be used so that i can compare the airlines based on each attributes. For example: if one attribute is "safety", i have end up in a result that "airline B is more safer comapred to other airlines (A,C,D,E,F)


title: Finding the most important product feature for a customer

text: Looking at a data of the format given below, I am asked to find the (top 3) feature(s) that affect a customer's purchase preference the most. In the dataset below, each customer row contains the break-down of purchases of a single product by various features of the product. `Fi` stands for the i'th product feature (e.g. colour, material) and `Catj` stands for the j'th value that a feature can take on (e.g. Purple, Woolen, etc). All features are categorical and they can take on different number of categories, e.g. `F1` below can take on 3 categories, `F2` can take on 2.                Customer  |  F1,Cat1  | F1,Cat2 | F1,Cat3 | F2,Cat1 | F2,Cat2 |...  | Fi, Catj|             1      |    100    |    800  |   100   |   1200  |   1200  |...  |   ...   |           2      |    300    |    300  |   300   |    100  |    900  |...  |   ...   |          3      |    250    |    250  |   250   |    200  |    100  |...  |   ...   |         ...    |    ...    |    ...  |   ...   |    ...  |    ...  |...  |   ...   |      I don't have a statistical background. However, I can visually interpret that Customer 1 is highly selective on `F1`, and majorly prefers `Cat2` over other categories of this feature. On the other hand, her purchases wrt Feature 2 show that she is not very selective about the binary `F2` feature. The situation is reversed for Customer 2. **Q1.** What is the most suitable formal method to carry out an analysis like this? (I thought this may be achieved by analysing the variances per feature.) **Q2.** Is it OK to assume a multinomial distribution for each categorical feature? If so, we get covariance matrices per each feature, and I am not sure how to interpret them. Any comments or pointers to relevant online sources will be greatly appreciated!


title: when can the Kalman filter return strange values?

text: Sorry if I'm asking a dumb question, but I'm not very good on statistics. I have tried to do an implementation of Kalman filter. I have double checked and the implementation seems correct with test values. But when I try it with real data (computer vision tracking system) I'm getting some `nan` or `inf` values.. So I thought maybe it could be wrong input data that are wrong. A colleague suggests that maybe Kalman filter fails if it is getting bad input like "if it get for several times a 0 value it will get infinite" So.. here is the question: is it possible that the Kalman filter gives as result an infinite number? and when?


title: Standard error of the quotient of two estimates (Wald estimators) using the delta method

text: I have two coefficients' estimates from a regression, each of which has an estimated standard error.   I would like to know the quotient of these two estimates -- that is, divide one of the estimates by the other. What would be the corresponding standard error? Would this be a candidate for the Delta Method? If so, how should the formula be applied? If this is not something that can be computed in a straightforward manner, is there a way to do this in Stata? In particular, I am interested in the Wald Estimator. In the Delta Method, there is a covariance term typically. Is it zero or non-zero in this case; that is, when can we assume that the two estimates are independent?


title: Non-linear (e.g. RBF kernel) SVM with SCAD penalties implementation

text: Is there one? I think there's a penalizedSVM package in R but it looks to use a linear kernel. Can't quite tell from the documentation. If it's linear, is there a R package that lets me calculate the kernel matrix so I can transform my data? I was thinking I could combine the two together to run a RBF SVM with SCAD penalties.


title: Which kind of statistical test should I try for a sample with bivariate data?

text: I'm trying to analyze data from a questionnaire, showing graphical results. I used dispersion charts for correlating bivariate data with continuous and discrete ranges, showing also a trend line on the distribution. What about correlating binary data samples and discrete samples? E.g., I have the following 2 series:               X = (1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0)      Y = (4, 4, 2, 4, 3, 4, 1, 4, 3, 3, 3, 3, 3, 3)      The first variable may have only values 0 or 1, the second one may have values in a range [1,4]. A friend suggested that t-test, Chi-squared and Fisher's test could be interesting. But I can't imagine how they can be used in a graphical way.


title: Calculating p-values and pnorm() in R

text: I am trying to calculate the p-values of observations by comparing them to the normal distribution in R using `pnorm()`. I have constructed a random distribution as my background model on which I would like to test the significance of various tests. I know for example, my background normal distribution has a mean of 1 and a standard deviation of 3. Say I have one test that I would like to test the significance of `test1 <- 20`. To obtain the p-value of a specific observation with a value of 20, I can use `pnorm(20, mean=1, sd=3)`. But what if for the same test, I have 5 repeated observations (technical repeats of the same test) with the values:               test1.rep1 <- 20     test1.rep1 <- 25     test1.rep1 <- 15     test1.rep1 <- 25     test1.rep1 <- 15      These 5 numbers have a mean of 20 and a standard deviation of 5. I could simply combine all of the repeated observed by taking the mean and then comparing to the normal distribution, i.e., `pnorm(20, mean=1, sd=3)`. But in this case, I am leaving out information, namely the standard deviation of these 5 observations. Is there an alternative way to include both the mean and standard deviation of the 5 observations when calculating the p-value? * * * The alternative is to calculate 5 p-values for each of the 5 observations of `test1`.               pnorm(20, mean=1, sd=3)     pnorm(25, mean=1, sd=3)     pnorm(15, mean=1, sd=3)     pnorm(25, mean=1, sd=3)     pnorm(15, mean=1, sd=3)      But then I have to find a way to combine these p-values, to end with one final p-value to look for significance. I ultimately want to know if `test1` is significant. I've looked into the Fisher method and the Stouffers method, but now I am thinking it may be better to just combine the values up front, rather than combining p-values.


title: Pre-process classification data according to ‘amount of evidence’

text: I have an area which is divided into polygons of different sizes. Each polygon has the same associated features/predictors and I know whether something occurs within the polygon or not. If something occurs in the polygon I have a number between 1 ... n which indicates the ‘amount of evidence’. I am trying to use logistic regression (or GAMs later on) to model the occurrence probability given a polygon’s features $(X_1, ..., X_n)$. At the moment I do not adjust for polygon size and ‘amount of evidence’ so I have just a vector $(X_1, ..., X_n)$ plus a zero or one value for the response. I am just wondering whether it may be useful to adjust the data (somehow?) to account of the different sizes of the polygons and the amount of evidence. Perhaps one should repeat the data rows more often depending on the size and the amount of evidence?


title: Stock Returns Covariance Prediction - Number of Principal Components

text: I am working on the following problem. Given N days of stock returns, I compute the covariance matrix for stocks. I then use Probabilistic PCA to "shrink" the covariance matrix. I am trying different number of principal components: 1, 2, etc. I then compute the log likelihood of observing the next-day's returns under the multivariate normal distribution using the "shrinked" covariance matrix (inverse of it to be precise). Just to make sure, the matrix is non-singular, as in PPCA the "remaining" principal components are not zeroed out, but instead are assumed to be uncorrelated. I run this process over some period of time, to compute the total log- likelihood of the next-day's returns prediction under this method (every day re-estimating the covariance matrix, computing next-day's returns log- likelihood, and adding it up). Just to make sure, the computed log-likelihood is always out-of-sample. Now, I would expect that the highest log-likelihood would be for a small number of principal components. However, what I find is that the highest log- likelihood is in fact for the model with 1 principal component, gradually decreasing as more principal components are used. Now, this could mean two things: 1) The best model is the one using 1 principal component (aka CAPM) 2) I need to adjust the computed log-likelihood to make it properly Bayesian. I know that if I would be fitting the model in-sample, then I would need to actually estimate the posterior for the model parameters (including the number of principal components), and then integrate it out. However, because I am calculating log-likelihood out-of-sample, it's in fact something like cross- validation. And so I don't really see the problem with this approach. Are there any problems I might have missed?


title: Confusion about the kind of offset-value for this non homogenous poisson-process

text: Suppose that I observe $M$-times this kind of table:                    t X      1:  1 0      2:  2 0      3:  3 0      4:  4 0      5:  5 0      6:  6 1      7:  7 0      8:  8 0      9:  9 0     10: 10 0     11: 11 0     12: 12 0     13: 13 0     14: 14 0     15: 15 0     16: 16 5     17: 17 1     18: 18 0     19: 19 0     20: 20 0     .   . .     .   . .     T-1:T-1 0     T:  T 0      So, what I get is like: $x_1 = 0, x_2 = 0, x_3 = 0, x_4 = 0, x_5 = 0, x_6 = 1....$, i.e., $X_t \in \\{0,1,2,3,4,5,...\\}$ where $t = 0,1,2,...,T-1,T$. As we can see, I observe jumps greater then one for particular $t$'s (like $t = 16, \rightarrow x_{16} = 5$). **I want to investigate if more counts occur if $t$ moves towards $T$ and what the impacts of different covariates, captured in $\boldsymbol{X}$, are.** To do that, I could assume that $N(t) = \sum_{i=1}^t X_t$ is a poisson- process, i.e., (1) $P(N(t) = k) = \frac{\lambda(t)t^{k}}{k!}\text{exp}\\{-\lambda(t)t\\}$ The regression model will be: (2) $\text{log}\\{E(N(t)|t)\\} = \boldsymbol{X\beta + \epsilon}$ $\Leftrightarrow \text{log}\\{\lambda(t)t\\} = \boldsymbol{X\beta + \epsilon}$ $\Leftrightarrow \text{log}\\{\lambda(t)\\} = -\text{log}\\{t\\} + \boldsymbol{X\beta + \epsilon}$ In this case, I would model $-\text{log}\\{t\\}$ as the _offset-value_ and take $N(t)$ as the dependend variable - which is just the cumsum of $X_i$, where $i$ goes from $1$ to $T$. **Question 1:** _I think here I'll get a problem because, as far as I remember, for particular $t$'s, the process is not supposed to jump by more then one?_ Another way I could think of, is, that I see my data as a big cross-table. If I assume that each cell count is a realization of an independent poisson- process and the total count $n$ is fixed, then I'm able to say that I observe a set of counts $n_1,...,n_T$ which are each poissonian with rate $\lambda_t$. If these assumptions hold, it can be shown that the cell counts are Multinomial$(\frac{\lambda_t}{n_1},...,\frac{\lambda_t}{n_T},n)$. Hence, if I see the observations as a multinomial respones, for the cell counts $X_t$, I'll get the model: (3) $\text{log}\\{E(X_t)\\}) = \boldsymbol{X\beta + \epsilon}$ $\Leftrightarrow \text{log}\\{\lambda_t\\} = \boldsymbol{X\beta + \epsilon}$ $\Leftrightarrow \text{log}\\{np_t\\} = \boldsymbol{X\beta + \epsilon}$ $\Leftrightarrow \text{log}\\{p_t\\} = -\text{log}\\{n\\} + \boldsymbol{X\beta + \epsilon}$ **Question 2:** _Now, I would choose $-\text{log}\\{n\\}$ as my offset-value. If the sequence of $\lambda_t, \ \ t = 1,...,T$ is increasing in $t$, this would answer my question. Since my $n$ is actually fixed and with this approach I'm able to deal with the multiple jumps at some $t$ it seems that the model (3) would be the way to go?_ **Question 3:** _If I use a`glm()` what would my dependend variable be? I thought about defining a new r.v. $Y_t$ which is 1 if $X_t>0$ and zero otherwise. If $X_t>1$ then I would duplicate the corresponding row $x_t$-times. If I do that, in my optinion, things are not supposed to change since I look at a ratio as a dependend variable anyway?_ **Question 4:** _Another issue I've is, that in model (2) I'm not sure how (3) catches that moving towards $T$ the cumsum of $X_t$ inceases steadily, since (3) looks at each cell separately._


title: significant predictors become non-significant in multiple logistic regression

text: When entered into univariate logistic regression I get the following:               Predictor 1: B = 1.049, SE = .352, Exp(B) = 2.85, 95% CI = 1.43 - 5.69, p = .003           Constant: B = -.434, SE = .217, Exp(B) = .648, p = .046          Predictor 2: B = 1.379, SE = .386, Exp(B) = 3.97, 95% CI = 1.86 - 8.47, p < .001           Constant: B = -.447, SE = .205, Exp(B) = .639, p = .029      and when entered into multiple logistic regression:               Predictor 1: B =  .556, SE = .406, Exp(B) = 1.74, 95% CI = 0.79 - 3.86, p = .171        Predictor 2: B = 1.094, SE = .436, Exp(B) = 2.99, 95% CI = 1.27 - 7.02, p = .012           Constant: B = -.574, SE = .227, Exp(B) = .563, p = .012      Both predictors are dichotomous categorical. I have checked for multicollinearity. Not sure if I have given enough info, but I cannot understand why predictor 1 has gone from being significant to non-significant and why the odds ratios are so different in the multiple regression model. Can anyone help (with a basic explanation of what is going on - I really struggle to understand statistics!!)


title: How to compute a new variable that is defined by two categorical variables using SPSS?

text: I have a questionnaire with a part (a) and part (b). It looks like this:               (a) "Do you ever feel sad":           0="never", 1="sometime", 2="often", 3="nearly always"     (b) "How distressed are you by this":          0="not distressed", 1="a bit distressed", 2="quite distressed",          3="very distressed".      I want to create a new variable:               "Do you feel sad?"     Response: 0="no", but 1="Yes, but no distress", 3="Yes and its distressing".      How can I create this new variable using SPSS?


title: Testing samples against a distribution

text: Suppose I have some data sitting in front of me, which is supposed to be sampled from a specific distribution. How would I go about testing whether it _is_ from that distribution? I've seen several tests which appear to do this - Pearson's chi-squared test, the one-sample Kolmogorov-Smirnov test, and the Anderson-Darling test, to name a few. Which one do I want here? (And if the answer is "it depends", then what does it depend on?)


title: How to determine if there is a drop in revenue after a change in the system?

text: We recently released a new redesigned e-commerce website for a client. They are claiming to be seeing a 40% drop in revenue due to the changes. We have daily sales data going back a few years. We are trying to determine if a recent drop in revenue after a site relaunch can be attributed to the redesign or just due to normal variance. What tools/methods can we use to either prove or disprove that the most recent decrease in revenue is statistically significant or not? Thanks in advance.


title: Data whitening for improving regression

text: My Data: $X_i= \\{0.4;~7;~1,000;0;~0.8;~1;~0;~40;~0.7;~1;~0;~89,100\\},~Y_i=345$ The training set size is $\approx35, 000$. $Y$ is the dependent variable and the task is to estimate its value (regression).   1. I normalized the $X$ data by scaling them and then applying a sigmoid function. I did the same with $Y$. I tried using a neural network with 3 layers:      * 1st input(neurons=number of features),      * 2nd - hidden=20 neurons,      * 3rd - output=1 neuron;      * each layer has sigmoid function),      * learning rate=0.9,      * momentum=0.5,  I also tried using `SVM(kernel="radial")` and `glm(family="gaussian")` models, but nothing performed well.   2. I applied a whitening procedure (to remove correlations) to $X,Y$ together as one matrix (by chance). And as the result the StdErr was less than 1%.  **So, the question is:** in my test data I have only $X$ and usually whitening is applied only to $X$ (inputs). Is it possible to apply results from previous whitening to a new dataset (without known $Y$)? Why did I get so much improvement?


title: What would be an illustrative picture for linear mixed models?

text: Say that you are in the library of your department of statistics, and that you come across a book with the following picture in the front page. ![enter image description here](http://i.stack.imgur.com/VWisF.png) You will probably think that this is a book about linear regression things. What would be the picture that would make you think about linear mixed models?


title: Can I adapt a MCMC proposal using a parallel chain?

text: I am running two MCMC chains (say chain A and chain B) in parallel, using the Metropolis-Hastings algorithm with acceptance probability: $P(accept\ x_t) = \min\\{1, f(x_t)/f(x_{t-1})\\}$. I would like my proposal distribution to be adaptive, but I know that if I use old values of one chain to choose the proposal (i.e $p(x_t|x_{t-1}) = g(x_1, \dots, x_{t-1})$), I will lose the Markov property (see this post for example), because $p(x_t|x_{t-1}, \dots, x_{1}) \neq p(x_t|x_{t-1})$. So I'm running two chains in parallel (same target distribution), chain A uses a fixed (non adaptive) proposal while chain B is using an adaptive proposal which is a function of the past values of chain A (i.e. $p(x^B_t|x^B_{t-1}) = g(x^A_1, \dots, x^A_{t-1})$). To me it looks like the algorithm is correct, because $p(x^B_t|x^B_{t-1}, \dots, x^B_{1}) = p(x^B_t|x^B_{t-1})$, given that the proposal for B is not using past values of the same. My questions are: 1) Is chain B still Markovian? 2) Even though the two chains have the same target distribution, is using one chain to adapt the other one a good idea? Maybe not because in the initial phase the two chains might be in different area of the parameter space. Thanks!


title: Connection between Fisher metric and the relative entropy

text: Can someone **prove** the following connection between Fisher information metric and the relative entropy (or KL divergence) in a purely mathematical rigorous way? $$D( p(\cdot , a+da) \| p(.,a) ) =\frac{1}{2} g_{i,j} da^{i} da^{j} + (O( \|da\|^{3})$$ where $a=(a^1,\dots, a^n), da=(da^1,\dots,da^n)$, $$g_{i,j}=\int \partial_i (\log p(x;a)) \partial_j(\log p(x;a))~ p(x;a)~dx$$ and $g_{i,j} da^{i} da^{j}:=\sum_{i,j}g_{i,j} da^{i} da^{j}$ is the Einstein summation convention. I found the above in the nice blog of John Baez where Vasileios Anagnostopoulos says about that in the comments.


title: Permutation tests in R for correlations

text: What's the easiest way to run permutation tests for a bivariate Pearson correlation matrix in R? I have tried running MPT.Corr, Multivariate Permutation Test for Correlations by Urbano Blackford et al., but I am not sure this function is applicable to the problem. The data contains mainly Likert scale items and time measures. At least one nominal variable needs to be entered as a covariate - i.e. controlled for.


title: Logistic / T-test p-values

text: I have Score data: integer ranging 0-100 which behaves normally on a QQ plot. The score is assigned to a subject which is either Normal or Diseased and it differs between two. I tried 2 models: \- prediction of Normal or Diseased status based on Score \- t-test of means between 2 groups Both gave me low p-value for the Score variable (t-test gave me the mean difference between 2 groups as well). Which result is more relevant to the following question: Does the score differentiate between Normal and Diseased? Thanks!


title: To gauge motivation score?

text: In order to gauge the level of motivations of my respondents to connect in Second Life, I have proposed to them to answer (with a Likert scale - 7 points: from totally disagree to totally agree) to 16 statements (items) that I have categorized myself in four motivation categories. For example the fulfillment motivations   1. The fact that my avatar gains a higher status (in terms of money, material possessions, reputation, etc) is important to me. (HIGH STATUS)   2. The fact that Second Life allows me to undertake and make money is important to me (START BUSINESS)   3. The fact that Second Life allows me to gain valuable knowledge in the virtual world (scripting language, etc) and / or in the real world (to follow courses in Second Life, etc) is important to me. (TO ACQUIRE KNOWLEDGE)   4. The fact that Second Life allows me to create whatever I want is important to me. (CREATE)   5. The fact that Second Life allows me to be altruistic (helping new residents, …) is important to me. (BE ALTRUISTIC) I would like to compute the average score of each respondents the fulfillment motivation                          HIGH STATUS|START BUSINESS|TO ACQUIRE KNOWLEDGE|CREATE|BE ALTRUISTIC|      Resp n°29        4             6             6          7         6        5,8     Resp n°30        2             4             6          6         4        4,4     Resp n°31        5             7             4          1         5        4,4      In place of computing a simple arithmetic average I envisage a principal component analysis If I do a PCA for the 5 assesments of the fulfillment motivations : Principal components/correlation Number of obs = 373 Number of comp. = 4 Trace = 5 Rotation: (unrotated = principal) Rho = 1.0000               --------------------------------------------------------------------------        Component |   Eigenvalue      Difference         Proportion      Cumulative     -------------+------------------------------------------------------------            Comp1 |      2.72909            1.81017              0.5458            0.5458            Comp2 |      .918928           .121757               0.1838            0.7296            Comp3 |      .797171           .242364               0.1594           0.8890            Comp4 |      .554806           .554806               0.1110           1.0000            Comp5 |  4.44089e-16            .                       0.0000            1.0000     --------------------------------------------------------------------------      Principal components (eigenvectors)               --------------------------------------------------------------------         Variable    |   Comp1   Comp2  Comp3   Comp4 | Unexplained      -------------+----------------------------------------+-------------     Statut_Elevé    |   0.2544    0.8212   -0.4932    0.1330 |            0      Lancer_Bus      |  0.5549   -0.3110   -0.2713   -0.1475 |            0      Créer   |   0.4279    0.0046    0.4411    0.7889 |            0      Altruisme   |   0.3693    0.3637    0.6442   -0.5625 |            0      Acquérir_C      |   0.5549   -0.3110   -0.2713  -0.1475 |            0      --------------------------------------------------------------------      I was wondering if I could not sort my items in both components:                Resp n° 29      Comp1:     (0.2544*4)+(0,5549*6)+(0,4279*7)+(0,3693*6) + (0,5549*6) = 12,8875        Comp2:   (0,8212*4) + (-0,3110*6) + (0,0046*7) + (0,3637*6) + (-0.3110 *6) = 1,7672      And after computing the mean of both components : 7,32735 Is this approach appropriate? If not, what can I do better than simple average items to calculate a score of achievement motivation?


title: Estimating an unknown restricted Markov Chain from partial measurements

text: There is an Markov chain $M$ defined on states $1, ..., N$ with the special property that it only has transitions $p_i$ from $i$ to $i + 1$ , $q_{i + 1}$ from $i + 1$ to $i$, and $r_i = 1 - p_i - q_i$ from $i$ to $i$ (for $i \in [N - 1]$; $p_N = 0$ and $q_1 = 0$ since there is no successor/predecessor state). We know our chain respects this special structure, but we don't know the transition probabilities. Our data about the process is generated as follows:   1. Start in some known state $n_0 \in [N]$    2. For $k = 0, ... , {t - 1}$:      1. Go to state $n_k$ in $M$     2. Take $n_k$ many steps in $M$     3. Set $n_{k + 1}$ as the current state in $M$   3. Output $n_0,...,n_{t}$ In other words, we only have partial information about our walk through M. But from many of these $n_0,...,n_{t}$ I want to be able to infer the transition probabilities of $M$. How should I go about this? Is there a standard procedure for this? If so, is there an implementation in Matlab (or R)?


title: Good source of information about AMOVA

text: I'm looking for a good source of information about Analysis of Molecular Variance (AMOVA), a kind of statistical analysis used in population genetics. Especially how are the $\phi$-statistics calculated and what is the difference compare to F-statistics (fixation index) in measuring the level of differentiation among diffrent population divisions.


title: Precise meaning of and comparison between influential point, high leverage point, and outlier?

text: From Wikipedia > **Influential observations** are those observations that have a relatively > large effect on the regression model's predictions. From Wikipedia > **Leverage points** are those observations, if any, made at extreme or > outlying values of the independent variables such that the lack of > neighboring observations means that the fitted regression model will pass > close to that particular observation. Why is the following comparison from Wikipedia > Although an **influential point** will typically have **high leverage** , a > **high leverage point** is not necessarily **an influential point**.


title: Cross-validation when there are near-duplicates in the data

text: I use cross-validation, with random partitioning, to test the generalization power of a certain classifier (a short-sentences classifier). The problem is, my dataset contains many almost-identical samples, such as "I do not accept your offer", "I can not accept your offer", etc. So, when I partition the dataset randomly, there are many samples in the train partition that are very similar to samples in the test partition, and I get exceptionally high performance, which is actually over-fitting. Is there a better way to test the generalization power of my classifier?


title: How to explain poor classification performance of recall when using SVM?

text: I applied SVM to perform the classification against several data sets. It turns out that the performance metric of recall is pretty bad for one data set. It has recall around 50% while other data sets have recall around 80%. For this kind of scenario, what are the possible approaches that are available to improve the recall? Besides, why some data sets can have a poor performance in terms of recall? How to analyze this kind of problem?


title: Mann-Whitney for groups filled Likert scale

text: My stats knowledge is less than basic, so apologies for possible confusion. I have 2 independent groups (N=14 and N=14) who filled in a Likert-scale based standardized questionnaire, consisting of 48 items. The scale measures types of religious coping, so I have 12 questions for each one of the 4 coping styles. I would like to compare those two groups: to see how differently they score on each one of coping styles. My questions are:   1. Is it OK to use the Mann-Whitney test? I suppose yes.    2. If yes, how do I enter data in SPSS: every response (1 for "strongly disagree", 2 for "disagree", etc.) of 28 participants (which will make a total of 336 entries) for each of 4 coping styles, then code them according to group and do the test? Or should I sum up the scores (which does not feel right, as far as I use a test specifically to treat ordinal data).    3. What is the best graph to present the results?


title: Regression Estimation difficulties

text: My regression problem is properly formulated, but is encountering serious computational difficulties. Dependent: $Y$ = multinomial Independent: $X_1, \dots, X_{90}$ = linearly independent set of variables. (I verified the independence. Afterall, I defined these variables). Consider design matrix $X$, Hessian $H$, and gradient $G$. **Difficulties:** condition_number($H$) = $10^9$ Variance Inflation Factors (VIFs): all of them $ < 5.5$ except for one variable which has $VIF = 15.6$ eigenvalues$(H)$: ranges from $-10^{10}$ more-or-less-smoothly to $-17.3$ This causes parameter estimation to go whacky - any sort of Newton-Raphson approxmation encounters numerical problems when computing $ H^{-1} \cdot G$ Any suggestions or ideas?


title: Determining butterfly weight distribution

text: I'm a researcher studying a specific butterfly species. I've conducted an experiment where I measure the weight of this species at different locations on the earth. The data is discretized because it is recorded to the nearest microgram. In some locations I have hundreds of measurements in other locations I have only a few measurements. How can I build a discrete distribution for overall butterfly weight from all of these locations? My specific problem is that I'm having a difficulty with determining how the effect of variable sample size should be integrated into my calculations. How can this data be combined to generate an accurate distribution of butterfly weight for this species? Thanks for the help.


title: Variable analysis in multiple linear regression

text: I'm investigating how some weather variables (15) affect electricity demand in a specific area during the last 20 years. I was thinking to perform the following steps: 1\. Perform Multiple Linear Regression on each subset of selected variables 2\. Save t-statistics (p-values) for each run Then, I would to show the statistics (median, min, max, quantiles) of the t-statistics for each variable in order to give an idea about which is the most influencing. Finally, I would also show the relationship between each variable and the mean square error obtained with regressions using it. Do you think this approach makes sense?


title: Different means with barchart and descriptives in SPSS

text: I have a dataset (in SPSS 22) with different measurements of children, with the children in rows and with the measurements as variables in the columns. For example, the activity levels of the children were measured on a couple of days. I have the variables `activity day 1`, `activity day 2` and `activity day 3` for which the values are on a continuous scale. I also have the variable `gender`. I wanted to know the average activity level on day 1, for both boys and girls. When I run the `explore`-function for `activity_day1`, grouped by `gender`, the mean activity for boys and girls is what I expected it to be (`7.36` for boys, and `6.71` for girls).               EXAMINE VARIABLES=activity_day1 BY GENDER       /PLOT BOXPLOT HISTOGRAM       /COMPARE GROUPS       /STATISTICS DESCRIPTIVES       /CINTERVAL 95       /MISSING LISTWISE       /NOTOTAL.      Then I made a grouped bar chart for the activity of boys and girls on the different days.               GRAPH       /BAR(GROUPED)=MEAN(activity_day1) MEAN(activity_day2) MEAN(activity_day3)     BY GENDER       /MISSING=LISTWISE.      The height of the first bar should represent the mean activity-level on day 1 for the different genders, right? But the first bar (`activity_day1`) for the boys reached all the way to a mean of `23` and for the girls it barely reached `1`. The other mean values for the other days are different as well. They are not what I expected them to be and are totally different from the mean-values I got with the `explore` function. How is this possible when I used the same variables? What did I do wrong?


title: Correct interpretation of Lmer output

text: I have produced the following model:               >lmer(TotalPayoff~PgvnD*Type+Type*Asym+PgvnD*Asym+Game*Type+Game*PgvnD+Game*Asym+                        (1|Subject)+(1|Pairing),REML=FALSE,data=table1)->m1          PgvnD=A percentage (numeric)     Asym= a factor 0 or 1     Type=a factor 1 or 2     Game= a factor 1 or 2      from this model the terms `Type`, `Game` and `PgvnD:Asym` were shown to be significant by removal from the model. `PgvnD` and `Asym` on there own were not significant but were left in the model because the interaction between them was. The summary of this model is as follows;               > m7     Linear mixed model fit by maximum likelihood      Formula: TotalPayoff ~ Type + PgvnD * Asym + Game + (1 | Subject) + (1 |Pairing)         Data: table1        AIC  BIC logLik deviance REMLdev      1014 1038 -497.8    995.6   964.4     Random effects:      Groups   Name        Variance Std.Dev.      Subject  (Intercept)   0.000   0.0000       Pairing  (Intercept) 716.101  26.7601       Residual              89.364   9.4533      Number of obs: 113, groups: Subject, 73; Pairing, 61          Fixed effects:                 Estimate Std. Error t value     (Intercept)   81.727      6.332  12.907     Type2          7.926      2.852   2.779     PgvnD         -8.466      7.554  -1.121     Asym1        -12.167      7.583  -1.604     Game2         15.374      7.147   2.151     PgvnD:Asym1   26.618      9.710   2.741          Correlation of Fixed Effects:                 (Intr) Type2  PgvnD  Asym1  Game2      Type2       -0.188                                 PgvnD       -0.218 -0.038                          Asym1       -0.620  0.081  0.189                   Game2       -0.483  0.009 -0.010 -0.015            PgvnD:Asym1  0.233 -0.267 -0.766 -0.328 -0.011      Am I interpreting these results correctly?   * `TotalPayoff` is higher when `Type=1` than in `Type=2`, it is also higher when `game=2` than when `game=1`.    * Also `TotalPayoff` increases significantly with `PgvnD` if `Asym=1` but not if `ASym=0` (indicated by significant interaction term but non-significant single terms). Also I notice that the `Subject` random effect has SD and variance of 0. Can this then be removed from the model? What does this really mean?


title: Change and outliers detection by means ARIMA (Tsay procedure)

text: I have a question about automatic on-line outliers and change point detection in time series data. Now I read the paper: "Outliers, Level Shifts, and Variance Changes in Time Series" Ruey S.Tsay It's written in the abstract: "..The methods employed are extremely simple yet useful.." But for me it's not very easy to understand these methods and especially to implement. If someone can give me a hint how to do that in R/Matlab etc I will be highly appreciated.


title: Probability distribution to simulate number of users

text: I need to simulate with probability distribution number of users who are listening radio during the day. Given data : max number of users is 1000, day has 24 hours (from 0 to 24), the highest number of users should be from 11 to 15. Which probability distribution will be the best and which parameters of distribution should I use?


title: Is there a function to check a set of events/random variables are mutually independent?

text: Say, 5 random variables A,B,C,D,E which can take values either 0 or 1. I have the theoretical population chances for the 32 cases of P(A,B,C,D,E). Obviously they add up to 1. How to check if A,B,C,D,E are mutually independent.                0 0 0 1 0   ->   0.039026946      0 0 0 0 0   ->   0.221152694      0 0 0 0 1   ->   0.000000000      0 0 0 1 1   ->   0.000000000      0 0 1 0 0   ->   0.081381771      0 0 1 0 1   ->   0.000000000      0 0 1 1 0   ->   0.014361489      0 0 1 1 1   ->   0.000000000      0 1 0 0 0   ->   0.126588070      0 1 0 0 1   ->   0.000000000      0 1 0 1 0   ->   0.022339071      0 1 0 1 1   ->   0.000000000      0 1 1 0 0   ->   0.046583025      0 1 1 0 1   ->   0.000000000      0 1 1 1 0   ->   0.008220534      0 1 1 1 1   ->   0.000000000      1 0 0 0 0   ->   0.174007266      1 0 0 0 1   ->   0.000000000      1 0 0 1 0   ->   0.030707165      1 0 0 1 1   ->   0.000000000      1 0 1 0 0   ->   0.064032769      1 0 1 0 1   ->   0.000000000      1 0 1 1 0   ->   0.011299900      1 0 1 1 1   ->   0.000000000      1 1 0 0 0   ->   0.099601970      1 1 0 0 1   ->   0.000000000      1 1 0 1 0   ->   0.017576818      1 1 0 1 1   ->   0.000000000      1 1 1 0 0   ->   0.036652435      1 1 1 0 1   ->   0.000000000      1 1 1 1 0   ->   0.006468077      1 1 1 1 1   ->   0.000000000


title: Local sensitivity around the maximum

text: After having inferred a statistical model for my physical system $Y=f(X_1, \dots X_n)$, I want to estimate the **local sensitivity** of each variable $X_i$ (one by one) **around the point of maximum response**. All the methods I know are based on the analysis of the first derivatives along each variables. Being located at a maximum, all of these derivatives are flattened, these approaches are irrelevant. When I plot $X_i \to f(X_1, \dots, X_n)$, with all $X_{j \neq i}$ fixed to $\arg \max f$, I can visualize the shape of the response, usually asymmetric bell curves, but I don't know how to quantify the sensitivity. Bonus points for an approach including sensitivity of joint variables :)


title: MS business intelligence development studio CART analysis

text: I'm doing some clinical database research and in an effort to lessen the burden on our statistical staff, I started to look for different software solutions to get the analyses I need, and that's where BIDS (business intelligence development studio). BIDS allows queries to be run against a SQL Server and store the results of those queries in a table or view. That table or view is then consumed by BIDS and logistic and linear regressions as while as CART analyses can be done on them.   BIDS Version The x axis that is cut off on the lift chart is 'overall population %' ![enter image description here](http://i.stack.imgur.com/LbpXf.jpg) A mining accuracy chart of the CART ![enter image description here](http://i.stack.imgur.com/Vm5Te.jpg) I'm not quite sure how this works in other stats packages, but in BIDS you define a certain percentage of your data set to train the model and the rest of the set is compared against that model and the lift chart shows the improvement in identifying the outcome you desire vs. a guess. I'm only vaguely familiar with CART analyses in the first place, and don't know the first thing about R, but this same analysis was done in R with very similar results. The red portion of what looks like a health meter in a video game corresponds nearly identical to the analysis done in R. However, there are no p values in the BIDS version. Consider the node Max Total Poly Pharm >=7 and < 13\. BIDS shows me (not present in the picture) value cases probability not present 2133 91.89 Is there any way to ascertain a p value from that. And does R use any of it's cases as training data for the model?


title: Does $y'y= \hat y' y + e'e$ hold for the least square model?

text: I want to show that $y'y= \hat y' \hat y + e'e$ hold for the least square model. I found out that: $\hat y= X b$ with $b$ the least squares estimator of the coefficient vector and $e$ the residual vector. I really appreciate your hints!


title: t distribution method of moments

text: This is a further question to my original question, where I did not get an helpful answer (at leas not helpful for me) :Methods of moments for t distribution I want to fit a t distribution to my data (data can be found here: http://uploadeasy.net/upload/bo9j6.rar). The probability density function is given by: \begin{align*} f(l|\nu ,\mu ,\beta) = \frac{\Gamma (\frac{\nu+1}{2})}{\Gamma (\frac{\nu}{2}) \sqrt{\pi \nu} \beta} \left(1+\frac{1}{\nu}\left(\frac{l - \mu}{\beta}\right)^2 \right)^{-\frac{1+\nu}{2}} \end{align*} First of all I use ML with the following code (Fitting t distribtution to financial data):               # fit t distribution     library(MASS)          fitdistr(alvsloss, "t")          # or          # log-likelihood function     loglik <-function(par){     if(par[2]>0 & par[3]>0) return(-sum(log(dt((alvsloss-par[1])/par[2],df=par[3])/par[2])))     else return(Inf)     }          # optimisation step     optim(c(0,0.1,2.5),loglik)      I get the following output:                  m               s               df            -0.0004919768    0.0130128873    2.6340459185       ( 0.0003182568) ( 0.0003453702) ( 0.1620424078)      and               $par     [1] -0.0004451138  0.0129659465  2.6182237477      which is more or less the same, I guess differences are due to the precision of the numerical procedures. Now I want to do this with using the method of moments, I use the following part of a paper (http://ideas.repec.org/a/eme/jrfpps/v7y2006i3p292-300.html): mean: \begin{align*} \mu=E(l) \end{align*} variance \begin{align*} \sigma^2 = V(l)= E((l-\mu)^2)=\frac{\beta \nu}{\nu-2} , \nu>2 \end{align*} excess kurtosis \begin{align*} \kappa=\frac{6}{\nu-4} , \nu > 4 \end{align*} my first questions:   1. Why are they using the excess kurtosis and not the third moment the skewness?   2. What values do I have to insert, is the following correct?: mean: \begin{align*} \mu=E(l)=\bar{l} \end{align*} variance \begin{align*} \sigma^2 = V(l)= E((l-\mu)^2)=\frac{\beta \nu}{\nu-2} = \frac{1}{n}\sum_{i=1}^n (l_i-\bar{l})^2, \nu>2 \end{align*} excess kurtosis \begin{align*} \kappa=\frac{6}{\nu-4} = \frac{1}{n} \sum_{i=0}^n \left(\frac{l_i-\bar{l}}{s}\right)^4-3, \nu > 4 \end{align*} this gives: \begin{align*} \hat{\mu}_{MM}=\bar{l}\\\ \hat{\nu}_{MM} = \frac{6}{\left(\frac{1}{n} \sum_{i=1}^n \left(\frac{l_i-\bar{l}}{s}\right)^4-3\right)} + 4\\\ \hat{\beta}_{MM} = \left(\frac{1}{n}\sum_{i=1}^n (l_i-\bar{l})^2\right) * \frac{(\hat{\nu}-2)}{\hat{\nu}} \end{align*} so I am using the sample mean, sample variance and sample excess kurtosis. Is this correct? And my main question: The output of ML tells me (the column df), that $\nu$<4 ($\nu$ is the number of degrees of freedom, df), but in MM I need $\nu$ to be greater than 4 or? So what does this mean? Is MM not usable? Or does it not matter?


title: How to choose kernel functions for support vector regression

text: Are there any good resources regarding how to design kernels for regression problems, specifically time-series regression type of problem. I am finding the choice of a kernel for regression extremely unintuitive. For a classification problem with hinge loss, the idea is that the kernel should map similar examples closely while keep dis-similar examples apart. What is the equivalent notion for the ILF loss used in SVR?


title: Data smoothing with derivative constraints

text: Which existing methods are capable of performing a smooth curve interpolation with constraints on derivative? What I need it to smooth my evenly spaced (with few missing points) data to get a smooth curve which will have positive first and second derivatives.


title: How do I combine multiple prior components and a likelihood?

text: Lets imagine I am comparing two groups of animals (treatment/control). There is previous data from cell cultures indicating the treatment should have a positive effect. This gives me "prior component 1". There are also two previous studies very similar to my own. One of them had an effect of 5 +/- 1 (prior component 2), the other of 1 +/- 2 (prior component 3). I feel the cell culture data is highly convincing, and that the prior component 3 is not such a reliable study. So I choose some weights of 3,1, and .5 for each and multiply. 1) To calculate the "overall prior" do I simply add these together as shown in the lower right panel? 2) Am I supposed to normalize these components before adding them? ![enter image description here](http://i.stack.imgur.com/U45TF.png) I then calculate a likelihood function for my current data as shown in the upper panel. ![enter image description here](http://i.stack.imgur.com/xSnyg.png) 3) How do I combine this information with the prior information shown in the first figure? For the lower panel I simply multiplied overall prior*likelihood. 4) I then want to make a decision based on this outcome. If I believe the effect is between -1 and 1 then I will stop studying the drug. If the effect is < -1 then I would perform new study A, if the effect is > 1 I will perform new study B. 5) Obviously there are a number of ways of choosing a decision (% density between -1 and 1, etc) Is there a best choice? 6) I feel I am doing something incorrectly, but maybe not. Is there a name for what I am trying to accomplish? **Edit** : If it helps I am trying to use the framework proposed by Richard Royall: 1) The likelihood function tells me "how to interpret this body of observations as evidence" 2) The likelihood function + priors tells be "what I should believe" 3) The likelihood function + priors +cost/benefit determines ""what I should do". Royall R (1997) Statistical evidence: a likelihood paradigm (Chapman & Hall/CRC) While the priors used here are subjective/nebulous they are built out of simple building blocks (uniform and normal distributions) that mathematically unsophisticated researchers can understand quickly. I think they convey my thought processes as a researcher well. Others may of course know of different background information. They should be able to build their own "compound prior" which may lead to a different decision than mine, but we should always agree on the likelihood function. This approach (if implemented correctly, which I am not sure I am doing here), appears to me to model the actual thought processes of researchers and thus be suitable for scientific inference. The steps map to the common sections found in scientific papers. The priors are the introduction, the likelihood is the results, and the posterior probability is the discussion. **R code:**               #Generate Priors     x<-seq(-10,10,by=.1)     y1<-dunif(seq(0,10,by=.1), min=-10, max=10)     y1<-c(rep(0,length(x)-length(y1)),y1)     y2<-dnorm(x, mean=5, sd=1)     y3<-dnorm(x, mean=1, sd=2)          #Weights for Priors     wt1<-3     wt2<-1     wt3<-.5          #Final Priors     y1<-y1*wt1     y2<-y2*wt2     y3<-y3*wt3          #Sum to get overall Prior     y<-y1+y2+y3          #Likelihood function for "current data"     lik<-10*dnorm(x, mean=1, sd=1)          #Updated Posterior Probability?     prob<-lik*y               par(mfrow=c(2,2))     plot(x,y1, ylim=c(0,1), type="l", lwd=4,           ylab="Density", xlab="Effect", main="Prior Component 1")     plot(x,y2, ylim=c(0,1), type="l", lwd=4,           ylab="Density", xlab="Effect", main="Prior Component 2")     plot(x,y3, ylim=c(0,1), type="l", lwd=4,           ylab="Density", xlab="Effect", main="Prior Component 3")     plot(x,y, ylim=c(0,1), type="l", lwd=4,           ylab="Density", xlab="Effect", main="Overall Prior")                         dev.new()     par(mfrow=c(2,1))     plot(x,lik, type="l", lwd=4, col="Red",          ylab="Likelihood", xlab="Effect", main="Likelihood")     plot(x,prob, type="l", lwd=4, col="Blue",          ylab="Probability", xlab="Effect", main="Posterior Probability?")     abline(v=c(-1,1), lty=2, lwd=3)


title: After centering, interactions' effects got reversed. It was not fixed by reversing the dependent variable. What should I do?

text: In a binary logistic regression, I centered all the independent variables (binaries, continuous, dependent, independent, all [I have 5 IVs: three binary and two ordinal]). Multicollineariy was eliminated. But all the effects got reversed. I could manually fix the report (by switching every + to a minus and vice versa in beta values) and in odds ratios, by dividing 1 by each odds ratio. But I didn't know if Wald values remained the same? So I reversed the zeros and ones of the dependent variable in order to simultaneously reverse again all the once-reversed betas (and fix them). Ok IV betas are now fixed. However, their interactions remain in the opposite direction! It seems that the interactions' directions depend on the IVs, and have apparently nothing to do with the dependent variable. I think I should temporarily reverse all the signs of the centered IVs so the beta signs for the interactions (and their odds ratios) get what should be. But in that case, the betas for the main variables would be again reversed! But I don't know if it is valid? Can I rely on such a model, in which the signs of betas are correct for the interactions, but reversed for main variables, Or a model in which main variables are correct but interactions are reversed? Perhaps I should reverse centered IVs together with the DV, in the hopes that betas get fixed for both IVs and their interactions. Do you have any solutions in mind? I would appreciate.


title: What statistical method used to compare mean and inspect developmental changes?

text: I use a multiple mediator model in describing my data, and age is the independent variable (IV). I want to see if there exist developmental changes at the level of my dependent variable (DV) that are related to variations of mediators 1 and 2. If I consider z-scores to standardize data between two age groups, then compare mean, some would say (on this website) that this type of method will ignore possible developmental changes in my data. Then, how to show that there exist developmental changes (either increasing, decrasing, constant)?


title: Cauchy M estimator of regression in R

text: Was wondering if anyone knows of an R package to estimate the Cauchy-M estimator of regression (see for example the end of this section, but with simultaneous estimation of the scale parameter as in section 2 of (1)). > (1) Mizera, I. Müller, C. H. (2002). Breakdown points of Cauchy regression- > scale estimators, Statistics & Probability Letters, Volume 57, Issue 1, > Pages 79-89.


title: Distance measure for multi-categorical responses

text: I have a data set of categorical data where each question can have more than one answer. This is a toy example: question one: what did you eat today?               subject 1 :  Potatoes, Apples          subject 2 : Apples          subject 3 : Honey, Potatoes, Apples      question two: which TV shows do you watch?                 subject 1 :  House, Dexter            subject 2 : TWD, Dexter            subject 3 : The news      I need to get a similarity measure between, for example, distfun(Potatoes and Apples, Honey and Apples). I dont think that ordinary multivariate categorical analysis measures deal with this kind of responses. Can someone share some light over which measures use by referring me to some explanatory paper that he knows please? Thanks in advance EDIT: 1.- The order of the responses is not important. 2.- I am not just looking for a way to solve the issue ASAP by using the a modified version of the "overlap" distance that is used for multivariate categorical data but to find the framework that deals with this kind of stuff. Thanks for your responses.


title: SPSS GENLINMIXED for crossed groups

text: thank you for any help you might offer. I have data in which the groupings are crossed, not nested. Each player responds to the same set of items. In MIXED, I give a separate intercept to each player and to each item -- but not to each combination of these -- using two separate command lines /RANDOM=INTERCEPT | SUBJECT(item) COVTYPE(ID) /RANDOM=INTERCEPT | SUBJECT(player) COVTYPE(ID) Now I want to do that in GENLINMIXED, and I can't see how. The DATA STRUCTURE option seems to require that the data is nested. As an experiment, I just ignored this and wrote the comparable two separate lines as above. But results not the same. Is it possible to do this in GENLINMIXED? Any reason why it would be possible in MIXED but not in GENLINMIXED? Thank you very much for any help at all!!


title: Statistical tools to be used for small data sets

text: I am a newbie to statistics. I am trying to model admission criteria for B-School based on one past years data. It has roughly 370 samples. I am trying to set admission criteria taking their placement performance (Placed or Not placed). I have their pre b-school performance as covariates. Roughly 85% of the people are placed. When I did binary logistic regression, most of the data is clustered around 0.6 to 0.8 with no clear separation to set a valid cut-off point. Hence I am getting very low accuracy for not place students (15%). I got one predictor (high school grades) as significant. What are other statistical techiques I can use for this small sample data.


title: TF-IDF for text classification by taking into account the document class

text: I am looking for a TF-IDF weighting for text classification (not document ranking/retrieval) which takes also into account the document class. For example let's use the typical spam/not spam classification example and suppose that we have a folder with spam documents and another folder with not spam documents. How could the TF-IDF take into account also the class of the document? I tried the most common TF-IDF formula and my accuracy with a simple SVM reduced to 20% from originally 63% (without feature selection with TF-IDF threshold).


title: Covariance of random vector multiplied with a random matrix

text: For a random vector $x$ multiplied by a non-random matrix $A$, $y=Ax$ the covariance matrix of $y$ is given by $\Sigma_y = E[Ax (Ax)^T] = E[Ax x^T A^T] = A E[x x^T ]A^T = A \Sigma_x A^T$, where $\Sigma_x$ is the covariance matrix of $x$. How can I calculate $\Sigma_y$ if $A$ is a _random_ matrix (uncorrelated with $x$) and $E[A]$ is known? The approach above does not seem to be directly applicable.


title: Intuition behind the t-distributions density function

text: I'm studying about Student's t-distribution and I started to wonder, how would one derive the t-distributions density function (from wikipedia, http://en.wikipedia.org/wiki/Student%27s_t-distribution): $$f(t) = \frac{\Gamma(\frac{v+1}{2})}{\sqrt{v\pi}\:\Gamma(\frac{v}{2})}\left(1+\frac{t^2}{v} \right)^{-\frac{v+1}{2}}$$ where $v$ is the degrees of freedom and $\Gamma$ is the gamma function. What is the intuition of this function? I mean, If I look at the binomial distribution's probability mass function, it makes sense to me. But t-distributions density function makes no sense at all to me...it is not intuitive at all at first sight. Or is the intuition just that it has a bell- shaped curve and it serves our needs? Thnx for any help :)


title: Can I learn SAS online?

text: I'll complete my Master's degree in a year and I was thinking of learning SAS as I've heard it gives you an edge over other candidates while being recruited. Could you tell me if I can learn it online and learn it well?


title: Is it possible to cluster standard errors on time variables?

text: Somewhat weird question: But is it theoretically possible to have panel data (many i, many t) where standard errors are clustered on t?


title: GLM with data piled up at zero

text: I am trying to run a model to estimate how well catastrophic illnesses such as TB, AIDS etc affect spending on hospitalization. Now I have "per hospitalization cost" as the dependent variable and various individual markers as independent variables, almost all of which are dummy such as gender, head of household status, poverty status and of course a dummy for whether you have the illness (plus age and age squared) and a bunch of interaction terms. As is to be expected, there is a significant--and I mean a lot, of data piled up at zero (i.e., no expenditure on hosptalization in the 12 month reference period). What would be the best way to deal with a model such as this? As of now I decided to convert the cost into ln(1+cost) so as to include all observations and then run a GLM model. Am I on the right track? Sorry for not posting about the cross-posting. I wasn't aware that was required. A small clarification- I am not treating cost as a count variable- I transform the cost variable as ln(1+cost) so that the observations where cost=0 do not drop out.


title: Posterior Marginal in Forward Backward

text: In computing Forward Backward Algorithm[http://en.wikipedia.org/wiki/Forward%E2%80%93backward_algorithm], it seems they are calculating posterior. I knew Forward Backward is an algorithm of Expectation-Maximization(EM) family and generally iterates over Expectation and Maximization steps, to converge. Is posterior getting same result, or is it a standard procedure to replace iteration or I am missing something to understand? If any one of the esteemed members may kindly guide me. Thanking you in Advance, Regards, Subhabrata Banerjee.