title: Gradient boosting algorithm (steps) question

text: So, far I have read following regarding boosting:   * Boosting is an ensemble technique.    * Train learner sequentially, where early learners fit simple models to the data.   * Analyze data for errors, that is, focus on records that have high errors and try to get them right.    * In the end, all the classifiers are given weights and combined to give a final prediction.  Gradient Boosting (Boosted tree) for regression: Repeatedly follow these steps: (1) Learn a regression classifier (2) Compute the error residual ( obs y - predicted y) per data point (3) Learn a new model to try to predict the error residual. So, the step I don't understand is, how do you combine step 1 and 3 ( and repeat the whole process again?) I have read that you add the two classifiers together, ( one from step 1 and second from step 3 and you repeat the process). So, my question is, how does this combining step 1 and step 3 take place and go to the next stage? Explanation with a few data points ( numerical example) to visualize the algorithm will be very helpful.


title: MADE and MSE pros and cons

text: When assessing the performance of an estimator, in which scenarios should one prefer the use of the **Mean Absolute Deviation Error** (MADE) over the **Mean Squared Error** (MSE) and vice versa? Edit / Clarification: Assume that we have i.i.d data $(X_{1},Y_{1}),\cdots,(X_{n},Y_{n})$. Let $(X,Y)$ denote a generic member of the sample who's conditional mean is denoted by $m(x)=E(Y|X=x)$. The performance of an estimator $\hat{m}(x)$ of $m(x)$ is often assessed either by: $$MSE(x)=E\left[\left\{ \hat{m}(x)-m(x)\right\} ^{2}|\boldsymbol{X}\right]$$ or by $$MADE(x)=E\left[\left|\hat{m}(x)-m(x)\right||\boldsymbol{X}\right]$$ where $\mathbf{X}=(X_{1},\cdots,X_{n}).$


title: Why is Pearson's correlation coefficient defined the way it is?

text: $$ r = \frac{{\rm Cov}(X,Y)}{ \sigma_{X} \sigma_{Y}} $$ I do not understand this equation at all. Where does it come from? From my personal understanding ${\rm Cov}(X,Y)$ comes from that fact that $X$ and $Y$ are dependent random variables, that is, $E[XY]$ is not the same as $E[X]E[Y]$. Is this analogous to saying that $P(A \cap B) = P(A)P(B|A)$ if $A$ and $B$ are not independent? I'm just confused as to why we want the ratio of $E[XY]-E[X]E[Y]$ over the product of the standard deviations for $X$ and $Y$.


title: Why is my replication of Silver & Dunlap 1987 not working out?

text: I'm trying to replicate Silver & Dunlap (1987). I'm just comparing averaging correlations or averaging z transform correlations and back transforming. I seem to not be replicating the asymmetry in the bias they find (back transformed zs are not closer to the population value for me than rs). Any thoughts? Is it possible that 1987 computing power just didn't explore the space enough?               # Fisher's r2z     fr2z <- atanh     # and back     fz2r <- tanh          # a function that generates a matrix of two correlated variables     rcor <- function(n, m1, m2, var1, var2, corr12){         require(MASS)         Sigma <- c(var1, sqrt(var1*var2)*corr12, sqrt(var1*var2)*corr12, var2)         Sigma <- matrix(Sigma, 2, 2)         return( mvrnorm(n, c(m1,m2), Sigma, empirical=FALSE) )         }      With these function it's easy to look at a bunch of correlations (basically replicate silver and dunlap 1987) and see the difference between averaging correlations and averaging z-scores and back transforming. Here's just one.               r <- 0.9     Y <- replicate(20000, rcor(10, 0, 0, 1, 1, r))     rs <- apply(Y, 3, function(x) cor(x[,1], x[,2]))     mean(rs) - r     zs <- fr2z(rs)     fz2r( mean(zs) ) - r      Just looking at the sample size of 10 and correlations of 0.1, 0.5, and 0.9 these are the results.                    rho  r bias   z bias          0.1  -0.006   0.006          0.5  -0.024   0.021          0.9  -0.011   0.011      And these are derived from Table 1 of Silver & Dunlap.                    rho  r bias   z bias          0.1  -0.007   0.003          0.5  -0.025   0.001          0.9  -0.011  -0.007      These are quite different results. From my test I'm seeing that it's just a matter of direction of bias, not magnitude. But, in the published paper they're finding much less magnitude with z. I couldn't find a published non- replication.


title: Bootstrapping for standard scores

text: I would like to compare data among different sets where the values span the same range (0,1] but will have different means and unknown distributions. Given the unknown distributions, I thought of applying bootstrap to estimate distribution paramaters, but my problem comes at how to standardize across sets. The way, that a z-score would be used for a normal distribution. I have been reading resources on the bootstrap, from which (I think) I understand that any parameter theta-hat can be determined with the resampling methodology, but then this will only provide an estimate of uncertainty in that parameter, is that correct? That is, if it were mean, I would get a standard error of mean which I could then not use to test individual values in the set (only means). If I adjusted by sqrt(n) to estimate the sample standard deviation, that would still not work for a non-normal distribution. I'm trying to do this in R. Idea thus far (which I hope is not the only way) is:   * determine an empirical distribution function for each resample    * calculate ecdf() output for a range (eg. 0 -1 in tiny increments)    * average these results across all resamples    * determine new ecdf from averaged results that can be used to "score" test values  Another idea or suggestion would be much appreciated, as I'm unsure whether the above is even valid. EDITS in response to comments (sorry wasn't clearer first time out)   * there are multiples (not just a pair of groups, but a larger set) in which I want to test (actually ten).   * the "standardization" is meant to measure a score. It's a combinatoric analysis where there are individual items which have scores for something that works in common among them. If I take three at a time, there is a greater likelihood that I will find something in common than if I try to find something in common among nine of them. As an example, choosing three gives a mean "score" of 0.175, choosing nine, a mean of 0.077.   * goal is that when performing some analysis, say a test in group three that yields a score of 0.2, is this "higher" or further from some measure of central tendency than a score of 0.09 for a test in group nine where the distribution will be different   * data are simulated prediction scores for an interaction and scale is arbitrary   * another way, perhaps, to phrase the problem in the context of bootstrapping: is it possible (or valid) for theta-hat, the parameter of interest to be estimated, to be an empirical distribution function? That is, the parameter of interest to be the distribution of the original population itself?


title: confidence interval for coefficient of variation

text: Does anyone know of an R package that computes upper and lower bounds for a coefficient of variation? Preferably using one of the methods described under distribution here: http://en.wikipedia.org/wiki/Coefficient_of_variation#cite_note-13 I can't guarantee that the data is normal ahead of time although negative and small positive values for the sample mean do occur with negligible frequency. Thanks!


title: Sequential hypothesis testing in basic science

text: I'm a pharmacologist and, in my experience, almost all papers in basic biomedical research use Student's t-test (either to support inference or to conform to expectations...). A couple of years ago it came to my attention that Student's t-test is not the most efficient test that might be used: sequential tests offer much more power for any sample size, or a far smaller sample size on average for equivalent power. Sequential procedures of varying complexity are used in clinical research but I've _never_ seen one used in a basic biomedical research publication. I note that they are also absent from the introductory level statistics textbooks that are all that most basic scientists are likely to see. My question is three-fold:   1. Given the very substantial efficiency advantage of sequential tests, why are they not more widely used?    2. Is there a drawback associated with the used of sequential methods that would mean that their use by non-statisticians is to be discouraged?    3. Are statistics students taught about sequential testing procedures?


title: Overlapping probability distributions

text: i am working with a set of distributions. I have so far analyzed several probability distributions with respect to each other, that means i have made a t-test, to see how probable it is that two events happen at the same time according to those distributions. That means i have computed:               let x1 and x2 be probability distributions         prob(x1=Z, x2=Z), z be any value      thats what the t-test gave me. though i used something similar so the distributions can have different variances. I have now computed the overlap of both distributions with respect to each other. what i want to do now is look at the whole thing backwards. I want to look at a certain range and be able to give a probability that there is 1 or two of those occurrences. perhaps another way to look at it is using dice. I have two dice with differing probability distributions (probably gaussian for me, but lets assume one throws 6 more often, the other one 3 more often). The first way to look at it would be to tell whether both dice throw the same number. What i want to do is to tell if both dice throw a number in a certain range, for example 5 and 6.               let x1 and x2 be probability distributions         prob(x1=Z, x2=Z), z be a value a certain range       I think it would help if you could tell me what the name of the thing is that i am trying to do, so i can read up. But i would also be happy to have a solution that works out of the box. any pointers how i can compute this? cheers and thanks for advice edit: provided an example.


title: Large deviations proof question

text: Below is part of the proof of large deviations result. K is cumulant generating function. Can anyone explain how the last step follows? ![](http://yaroslavvb.com/upload/large-dev.png) This is page 157 of McCullagh's "Tensor Methods in Statistics"


title: Order values of $x_1$ to achieve desired correlation with $x_2$.

text: Let $D$ be a data set containing $n$ instances of two variables $x_1$ and $x_2$. Further let $x_1$, $x_2$ be independent according to some measure like linear correlation or mutual information etc. How do we find a reordering of the values of $x_2$ that results in some predetermined (achievable) correlation coefficient against $x_1$? I guess that a simple way to achieve (close to?) maximum correlation would be to order both in increasing order (independently) and then reorder row-wise into $x_1$'s original positions. Thanks


title: How do I estimate a repeated measures ANCOVA in R and SPSS?

text: I'm currently running a perception experiment:   * DV: error (this in degrees- how much an observer was away from the real answer)    * IV: the time bin (5 levels) in which the unique stimulus appeared   * IV2: True/False- whether the unique stimulus occurred before or after stimulus2.    * Covariate: Distance from the fixation point this stimulus appeared (continuous) So, I have decided to use repeated measures ANCOVA since, all the observers were exposed to all 5 levels of the IV multiple times. Currently I have written:               data.aov <-aov(error~(timebin*after*distance) +                     Error(subject/timebin*after*distance), d)     summary(data.aov)      **Is this the correct way to specify my repeated measures ANCOVA in R?** Also, I wanted to run this analysis in SPSS to check that I get the same results. However, SPSS doesn't like the format of my data. For repeated measures analysis, the 5 bins should be in 5 different columns, but in my data file, they are all together under the same heading `timebins`. **How can I run the repeated measures ANCOVA in SPSS if my data are in long format?** * * * Thanks everyone for great answers! @Marcus, just one thing- regarding regression. If I were to use regression, would I be looking at trends? I am actually more interested in comparing bin 1 and bin 5. This is why I was going to use a t-test, but now it seems like I am using within-subjects anova. Also, I've had a look at distance and (as you said) might not be a covariate, since it was randomised for every condition! Then the codes should look like this?               data.aov <-aov(error~(timebin*after) + Error(subject/timebin*after), d)     summary(data.aov)      But what is the difference between using `*` and using `+` ? Should I use `*` to get an interaction effect?


title: Appropriate significance test for peak level in a time series

text: I have collected biometric data on two stimuli (they are TV ads). The output is a second by second response, with values from $0$ to $1$. I want to know whether the peak response from one ad is higher than the peak response from the other. What is the appropriate test to perform to determine whether the difference is statistically different? I know that if I averaged the response over the entire ad I'd use a $t$-test to see if there was a difference, but I'm uncertain how to look at this for a specific point in time.


title: Difference between pairwise ANOVA and pairwise t-test?

text: Hellow I need to know what the difference between a pairwise t-test and a pairwise ANOVA is. Preferably explained in laymen terms.


title: Confusion related to calculation of partial derivative

text: I have this function $P = f(\alpha)$. $\alpha$ is a function $\alpha = f(\theta, x)=\theta x$. Now I have $P = \alpha(x_1+x_2 + ...x_n)$ Now I need to calculate the partial derivative of P wrt $\theta$. Then $\frac{\partial P}{\partial \theta} = \frac{\partial P}{\partial \alpha} * \frac{\partial \alpha}{\partial \theta} = (x_1+x_2+x_3)*\frac{\partial \alpha}{\partial \theta}$ Now what I am trying to do is use gradient descent to maximize P wrt $\theta$. So what I was expecting was the partial derivative wrt $\theta$ to be a single value. However, Here $\frac{\partial \alpha}{\partial \theta} = \frac{\partial \theta x}{\partial \theta}$ changes wrt to x. So I am a bit confused how to get a single value out of it?


title: Nonconstancy of error variance

text: If the residuals show that the non constancy of the error variance is clearly present, does it mean that your regression results are completely invalid?


title: Validating principal component analysis

text: I just wanted to do this small experiment to make sure I understand PCA correctly. My dataset contains 8 columns. The first two columns are randomly generated in excel => randbetween(4, 5) and the other 6 columns are also generated in the same way but the formula used is => randbetween(1,3) When I do PCA on this I am not getting good results. I expect that the result should indicate high eigenvalues for a factor that is a combination of first two columns and low on other colums. This is my code in R :               sensex.dat = read.csv('C:/Study/_SEM4/brand man/emperical/dice.csv', header = T)     attach(sensex.dat)     sensex.cov = cov(sensex.dat)     sensex.eigen = eigen(sensex.cov, symmetric = T)     sensex.eigen$values     sensex.eigen$vectors


title: What questions is your data team hoping to answer?

text: I'm interested in hearing about what statistics problems other data teams are trying to tackle. For me, I have data on energy usage in the United States, and we are trying to figure out ways to forecast energy demand as well as figure out ways to balance heavy loads on the power grid. The statistical challenges are in short term and long term time series forecasting. Regime changes such as weekdays -> weekends further complicate the problem, as well as external factors like whether variables. Spatial correlation is also a nontrivial issue. What type of data do you have, and what questions are you hoping to answer? What statistical challenges are there?


title: Local linear regression: number of grid points?

text: Local linear regression is a popular tool. How can one choose the number of grid points for which to estimate the unknown function on?


title: Proving the convergence of KDE algorithms when the samples are non-i.i.d

text: I am currently working on convergence proof for a new method for non- parametric importance sampling, and I need some help... My method uses an MCMC algorithm to generate a set of dependent $M$ samples $X_1 \dots X_M$ from a distribution $g$, and then reconstructs an approximation of this distribution $\hat{g}$ using a KDE algorithm. I want to show is that $\hat{g} \rightarrow g$ as $M \rightarrow \infty$. Foruntately, I was able to find a proof that is similar to the one that I need in a paper about non-parametric importance sampling by Ping Zhan (see the appendix in this paper here for more details ). However, in this paper, it seems as if the author assumes that the samples used in the KDE algorithm $X_1 \dots X_M$ are i.i.d - which does not hold in my case given that I used an MCMC algorithm to generate them. I am wondering if there is a way to tweak this existing convergence proof to account for the fact that we are using non-i.i.d samples within the KDE algorithm? Intuitively, I think that using correlated samples within a KDE algorithm will reduce its effectiveness - so maybe there is some kind of bounding factor? Otherwise, perhaps I exploit the fact that the samples come from an MCMC algorithm whose stationary distribution is $g$? If it helps, I can assume that the MCMC algorithm is just generic Metropolis-Hastings. I'm open to any ideas!


title: Validity of pseudo-panel data constructed from repeated cross sectional data as a panel data

text: I am looking at the repeated cross-sectional data from federal reserves, which has both panel data and repeated cross sectional data at different time- points,e.g. 2007-2009 is a panel while 2010 is a cross sectional data set and everything before that is repeated cross section as well until you get back to the 1983-1989 period which is also a panel. I want to use recent data-sets like 2001 - 2009, of which only the last two years will be true panel data. RCS data is considered to be inferior to true panel data in general in the sense that in the former case, the same individuals are not followed over time, thus making individual histories unobtainable to include in a model. However, Several authors such as Deaton (1985), Moffitt (1990,1993) showed that the RCS data can be used to estimate a few commonly used models such as the fixed effects model or the linear dynamic model. These methods are based on grouping “similar” individuals in cohorts and the ‘cohort-averages’ are treated as observations from a pseudo-panel. Note that, all the prior studies were conducted on repeated cross-sections without the panel part. Now, my first question is, 'Is there any known method to compare pseudo-panel data to a panel data?'. My idea is to fit a model to the synthetic panel data and estimate the parameters, then fit the same model to genuine panel data and compare the estimation accuracy. Does this sound correct? Of course, I want ideas about how much of it is doable. (Please note that I have limited ideas about how to manage a huge data-set like the ones available in fedres website.) Apologies if this was an irrelevant or ill-posed question.


title: Making new variable instead of correcting for temporal autocorrelation in a GLMM. Is it a valid alternative?

text: I am doing some forest disturbance research, in which the aim is to predict the probabilities of wind damage occurrence in forest stands of different site (altitude, slope steepness) and stand properties (age, timber volume). I am using a logistic regression (because of the binomial response variable: 0-1) with the R package "lme4". My data set looks like this (several thousand forest stands in the full data set):               stand_id    year    damage  altitude    slope   age volume  occurrence     123 2001    0   900 15  100 235 0     123 2002    0   900 15  101 242 0     123 2003    0   900 15  102 249 0     123 2004    3.6 900 15  103 256 1     123 2005    0   900 15  104 259.4   0     123 2006    2.1 900 15  105 266.4   1     123 2007    1.8 900 15  106 271.3   1     123 2008    0   900 15  107 276.5   0     123 2009    0   900 15  108 283.5   0     123 2010    0   900 15  109 290.5   0     124 2001    0   1100    10  80  172 0     124 2002    0   1100    10  81  181 0     124 2003    6.2 1100    10  82  190 1     124 2004    8.9 1100    10  83  192.8   1     124 2005    2.4 1100    10  84  192.9   1     124 2006    0   1100    10  85  199.5   0     124 2007    0   1100    10  86  208.5   0     124 2008    5.5 1100    10  87  217.5   1     124 2009    2.4 1100    10  88  221 1     124 2010    0   1100    10  89  227.6   0      . . . If damage is bigger than 0, the response variable (occurrence) has a value of 1, otherwise 0. The model would be: occurrence ~ altitude + slope + age + volume As you can see, there are measurements in all years for a 10-year period in all of the stands. Since we have more than one measurements for the same forest stand, a random effect (of forest stand) has to be added to the model, which will be the factor "stand_id". Regard to my understanding of Statistics, spatial autocorrelation is not a problem anymore, but temporal autocorrelation of the (residuals of the) response variable is still something needs to be dealt with. I have checked the autocorrelation function (acf) of it, and it is apparent in the first following year (acf value is ~ 0.1), then it dies off. Many experts advise to model the temporal autocorrelation structure and use it as a correction for the previously fitted (logistic regression) model. I find this a bit tedious and overcomplicated (also cannot do simply in lme4), therefore my question would be: Is it a valid alternative if I make a new predictor variable that is the damage of the previous year (because of the acf) and use it in the model, as well? To make it more clear, here is the new data set (last column is added):               stand_id year damage  altitude slope    age   volume occurrence dmg_prev     123 2001    0   900 15  100 235 0   0     123 2002    0   900 15  101 242 0   0     123 2003    0   900 15  102 249 0   0     123 2004    3.6 900 15  103 256 1   0     123 2005    0   900 15  104 259.4   0   3.6     123 2006    2.1 900 15  105 266.4   1   0     123 2007    1.8 900 15  106 271.3   1   2.1     123 2008    0   900 15  107 276.5   0   1.8     123 2009    0   900 15  108 283.5   0   0     123 2010    0   900 15  109 290.5   0   0     124 2001    0   1100    10  80  172 0   0     124 2002    0   1100    10  81  181 0   0     124 2003    6.2 1100    10  82  190 1   0     124 2004    8.9 1100    10  83  192.8   1   6.2     124 2005    2.4 1100    10  84  192.9   1   8.9     124 2006    0   1100    10  85  199.5   0   2.4     124 2007    0   1100    10  86  208.5   0   0     124 2008    5.5 1100    10  87  217.5   1   0     124 2009    2.4 1100    10  88  221 1   5.5     124 2010    0   1100    10  89  227.6   0   2.4      . . . Running the model with this new variable (dmg_prev), it also turns out to be significant, and if I am right, it carries the information on the temporal autocorrelation, as well. So no additional correction is needed anymore. Could someone verify this? Thank you very much for your help in advance!


title: Why sparse QR decomposition yield different results from dense QR decomposition in R?

text: A simple example is given below:               mat <- diag(1:6)     qr.R(qr(as(mat,"Matrix")))     qr.R(qr(dd))      the diagnal values are of opposite signs. There is a warning message when doing sparse QR decomposition:               Warning message:     In qr.R(qr(as(dd, "Matrix"))) :     qr.R(<sparse>) may differ from qr.R(<dense>) because of permutations      However, if I check `qr.R(qr(as(mat,"Matrix")))@q`, it is `0 1 2 3 4 5`, it seems that there is no permutations being done. Moreover, in Matlab, if I want to do a sparse matrix decomposition:               mat=sparse(1:6,1:6,1:6)     p = colamd(mat)     R = qr(mat(:,p))      this result(`R` and `p`) should be equal to the R output `@R` and `@q`, am I right?


title: Regression with a ordinal scale dependent variable in the context of panel data

text: In my panel data I observe a quite large amount of different individuals for a fixed time duration. The time axis is fixed for every individual, i.e., I observe 20 time points for every individual. The dependend variable is of kardinal scale but I choosed to transform it into ordinal scale. The reason for this transformation is due to the fact that I observe the case "-Inf" very often which can obviously not be included in any form of regression model. So I choose to split the axsis into intervals like (-Inf,-20000],(-20000,-10000],(-10000,0],(0,10000] and so on. Due to this transformation I actually observe a dependend variable of ordinal scale. My first question in this context is: In what way a ordinal scaled dependend variable with time and individual specific variation can be modelled. 2nd: If I acutally do a "ordinal regression" how do I deal with a continuous explanatory variable. I already checked how small the cells get and I'll receive a very large multidimensional table with a lot of zeros and here and there a number from 4 upwards. Kind regards


title: How to compute the generalization error (test error) in a multivariate regression?

text: How can I, after running multivariate linear regression, compute the Generalization Error (also known as Test Error)? Are there many practises to do that? EDIT: Can the following be considered as Generalization (or Test) Error?               y = data1; [M,N] = size(y);     X = [ones(M,N) data2];     y_train = y(M-30-483:M-30,:);     X_train = X(M-30-483:M-30,:);     y_test = y(M-29:M,:);     X_test = X(M-29:M,:);          [b,~,r,~,stats] = regress(y_train,X_train);     [M,N] = size(data2);     figure; plot(r);     axis([0,N+1,-Inf,Inf]);          y_estimated = X_test * b; y_expected = y_test;     generalization_error = y_expected - y_estimated;     figure; plot(generalization_error);


title: Including non-hypothesized interaction terms in moderation analysis

text: I am running a moderation analysis using hierarchical linear regression. I have two predictor variables and two moderating variables and I am interested in their interactions to predict an outcome. My analysis includes the four possible interaction terms, but I am only theoretically interested in two - is it ok to leave the other two interaction terms in analysis or should I only include the interaction terms that I specifically predicted?


title: Data mining techniques in R for advertising and sales data

text: I would like to apply one or more data mining techniques to a dataset, in order to see the effect advertising has on sales. I am working from this dataset. It has 36 consecutive entries of monthly data for both sales and advertising - it's very small. I exported the dataset to a ".csv". I deleted the date column, because I will use R's ts (time series object). The ".csv" now looks like this:               Advertising,Sales     12,15     20.5,16     21,18     ..., ..., ...     23.4,17     16.4,1      The example coded below works. However, I had to split the matrix into two lists, because of the `HoltWinters()` function. I would prefer to analyse Advertising and Sales together at the latter stages. What other data mining techniques may be more beneficial?               data <- read.csv("./advertising_sales.csv", header=TRUE)     data_ts <- ts(data, start = c(2011,1), frequency = 12)     print(data_ts) # to check data has been correctly added          > Jan 2012    13        17.3    21     + Feb 2012    14        25.3    29       + ...     + Nov 2013    35        23.4    17     + Dec 2013    36        16.4     1          plot(decompose(data_ts))     data_ts_ad <- data_ts[,1] #assign advertising as list, for HoltWinters     data_ts_sa <- data_ts[,2] # assign sales as list, for HoltWinters          #do HoltWinters for advertising     plot(HoltWinters(data_ts_ad))     data_ts_ad.hw <- HoltWinters(data_ts_ad)     predict(data_ts_ad.hw,n.ahead=9)          >           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug     + 2014 18.52852 25.47521 27.16683 36.41340 38.14678 33.04452 33.22488 32.12758           Sep     + 2014 32.58964          plot(data_ts_ad,xlim=c(2010,2014))     lines(predict(data_ts_ad.hw, n.ahead=24), col=2)          #do HoltWinters for sales     plot(HoltWinters(data_ts_sa))     data_ts_sa.hw <- HoltWinters(data_ts_sa)     predict(data_ts_sa.hw,n.ahead=9)          >          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug     + 2014 11.05723 23.27877 50.06859 57.22696 61.50669 26.35195 62.26159 70.83347           Sep     + 2014 23.18957          plot(data_ts_sa,xlim=c(2010,2014))     lines(predict(data_ts_sa.hw, n.ahead=24), col=2)      I recently came across a book called R and data mining: Examples and Case Studies by Yanchang Zhao. It has excellent worked examples and this is where I have found inspiration. However, I can't get my small brain to think which techniques can be applied to this dataset. I am new to R, so please try dumb-down your answers slightly. **EDIT:** Output of data_ts is given below.               dput(data_ts)          structure(c(12, 20.5, 21, 15.5, 15.3, 23.5, 24.5, 21.3, 23.5,      28, 24, 15.5, 17.3, 25.3, 25, 36.5, 36.5, 29.6, 30.5, 28, 26,      21.5, 19.7, 19, 16, 20.7, 26.5, 30.6, 32.3, 29.5, 28.3, 31.3,      32.2, 26.4, 23.4, 16.4, 15, 16, 18, 27, 21, 49, 21, 22, 28, 36,      40, 3, 21, 29, 62, 65, 46, 44, 33, 62, 22, 12, 24, 3, 5, 14,      36, 40, 49, 7, 52, 65, 17, 5, 17, 1), .Dim = c(36L, 2L), .Dimnames = list(         NULL, c("Advertising", "Sales")), .Tsp = c(2006, 2008.91666666667,      12), class = c("mts", "ts", "matrix"))


title: Adding predicted probabilites from logistic regression instead of using cut value

text: I am using a logistic regression model to predict a binary decision (purchase, don't purchase) based on several independent variables (income, age, education, etc.) for a population of individuals (customers). I have data for individuals from one or more previous time periods, and I want to predict behavior for different individuals in a future time period. Unfortunately, my experience is with explanation, not prediction. My real interest is in predicting aggregate behavior--for example, what are total predicted purchases by customers in a future time period based on their characteristics? I can see two ways of doing this. First, I could use the parameters from the logistic regression model to generate a probability [0-1] for each customer in the future time period, then use a cut value (0.5) to resolve those probabilities to either 0 or 1, then sum the 1s to generate an estimate of total purchases. Second, I could use the parameters from the logistic regression model to generate a probability [0-1] for each customer in the future time period (as before), then simply sum those probabilities to generate an estimate of total purchases (without using a cut value). The second approach (adding the probabilities) makes the most sense to me, but the reference material I have consulted so far frames the prediction task in terms of cut values and classification tables. Is the second approach conceptually flawed? If so, why? Thanks very much. ADDENDUM: With regard to the references I consulted, it was often suggested to use cut values and classification tables, with training and validation sets, to evaluate the real-world performance of a logit model. However, I would have thought that summing the probabilities would have been a better way to do that.


title: Query on DBN Implementation

text: I am Pradeep pursuing M.tech at S.I.T, Tumkur, India. I have stuck in a problem while implementing Deep Belief Networks for Isolated Word Recognition. The procedure i have adopted is described below, please do help me in analysing the stage where i went wrong. Back-propagation algorithm implementation : The Normal back-propagation algorithm that I have used for Isolated Word Recognition contains the following details of the feed forward architecture;   1. Input matrix -180×132 i.e., 180 inputs are fed to the FFNN and 33 such vectors belonging to each class. So there are 4 classes hence 33×4=132.    2. I have used a single hidden layer with 30 hidden nodes and the activation function is ‘sigmoid’ (1/1+exp(-x))   3. The output layer has 4 nodes so as to classify 4 different classes with the same “sigmoid” activation function.   4. The target labels for 1st class is [1;0;0;0] which is repeated 33 times, 2nd class has [0;1;0;0] repeating 33 times , 3rd class has [0;0;1;0] and 4th class has [0;0;0;1] each repeating 33 times.  The procedure that I followed for classification is through usage of Back- propagation learning Algorithm. Once the training is completed, I give a test input of dimension 180×1, correspondingly find the output by using the updated weights, the position corresponding to minimum mean square error between expected and obtained output is said to be the recognized word. The Classification result is approximately 92%. Deep Belief Network Implementation Now my intension is to perform the same task using Deep Belief Networks 1. Initially, the weight initialization is random but 1st pair of layers (input- hidden1) is treated as an RBM and the RBM update procedure is trough contrastive divergence procedure. 2. The procedure that I have followed is referred from Yoshua Bengio in “Learning Deep Architectures for AI, 2009”.   1. But the changes that I have made here is instead of using “sigm” at visible layer, I have used “normpdf” following Gaussian distribution with 0 mean and unit variance . This is because for real valued speech features, its necessary to use Gaussian distribution as provided in the literatures.   2. Once the first layer RBM is trained using the above algorithm, the activations at the hidden layer are fed as the input to the next layer (stacking RBM to form DBN) and is pre-trained in the same fashion as that followed by first layer. Now here are my further steps that I have implemented with lot of doubts.   3. I have added the output layer with 4 nodes so as to classify 4 different classes with “softmax activation function”   4. The cost function here is “Cross entropy” instead of “mean squared error” between the target (as illustrated in step 4 of back –prop algorithm) and the obtained output.   5. But from the results the error increases after every epoch (iteration) and becomes “NaN” (undetermined or undefined numerical results) after a few iterations. Here are my following queries about the above implementation   1. Is the procedure adopted is correct?   2. Logically I have used 2-hidden layers in the above described DBN implementation. Is that step correct? Because in the normal Back-prop algorithm I have used only one hidden layer.   3. Should I add immediately a output layer after pre-training 1st RBM? I kindly request you to give suggestions in correcting the implemented topics.


title: How is this version of the pearson correlation coefficient derived?

text: After taking my first semester stats class I started getting more and more into statistics. Currently I've come across this version of the Pearson Correlation Coefficient. Comparing this formula to the formula presented on Wikipedia shows they're inconsistent. Further, I've seen at least 3 different ways to calculate it. Here's the link: http://www.stat.wmich.edu/s216/book/node122.html and the formula therein: ![enter image description here](http://i.stack.imgur.com/DW9Np.jpg) How is this formula derived? Since it is dividing through by n I assume this particular version is for populations and also for discrete data. Am I correct?


title: Stats test which test of the NULL that a distribution is uniformally distributed

text: I'm looking for a statistical test which tells the probability that a given sample comes from a uniform distribution. Shapiro test wether a sample comes from a normal distribution. I'm looking a similar test which test wether a sample comes from a uniform distribution. Thank you


title: Correlation on several independent populations

text: I am performing a study in which i want to analyze the relationship between software project metrics. In each project, there is a population of software modules. For each module, I compute the values of two software metrics A and B. I want to compute rank correlation between the two metrics A and B, and I am not sure of the best way to do that:   * Should I perform a different correlation test on each project in the sample?    * Or is it possible to simply consider all the modules from all the projects together, and perform a single correlation test? By doing this I am afraid I would miss some information (e.g. A and B are correlated in some projects, but not in the others)


title: One step ahead forecast with SEASONAL data collected sequentially

text: In this post it was asked how to do one step ahead forecasts using Arima form the forecast package. Now I'm using an example with hourly seasonal data and would like to do something similar but the forecast will be 24 hours ahead. When I get new 24 hours of data I will add them up and produce another 24h forecast. Again without re-estimating the model. I'm not sure whether the following code is right since the fitted values are theoretically only for the next hour but not for the next 24 hours:               library(expsmooth)     data(utility)     n=length(utility)     y=ts(log(utility[1:(n-28*24)]),f=24)     new.data=ts(log(utility[(n-28*24+1):n]))          library(forecast)     model = auto.arima(y)     newfit <- Arima(new.data, model=model)     onestep.for <- fitted(newfit)


title: What impact does including more predictors have on the CI for R squared?

text: More precisely I mean " _All other things being equal_ , what impact does including more predictors ('IVs') have on the CI for R squared?" Say we'll get an R squared of .30 no matter whether we include 4 predictor variables or 7 predictor variables. Would one definitely produce a more precise CI than the other, or tend to produce a more precise CI than the other? Or is the number of predictor variables unrelated to the precision of the CI?


title: Mechanics of and intuition behind probabiliity

text: Suppose in a given experiment that an event $E$ will occur with probability $p$; if the experiment is staged only once, then we know from the previous statement that the probability that the event $E$ occurs is given by $p$. Now, say we stage the experiment $n$ times (where $n \to \infty$ and experiments are independent); intuitively (and mathematically, too, I guess), because the probability is constant across "stagings" of the experiment, I have deduced (perhaps erroneously) that the probability (before any staging has taken place) of the event $E$ occurring once across $n$ stagings is also given by $p$. However, this seems to conflict with the frequentist definition of probability, paraphrased by the following: Given $p = \frac{x}{y}$, as $n \to \infty$, it should be expected that, on average, event $E$ will occur $x$ times per every $y$ iterations of the experiment. Consider an example: A man with a deck of cards (52 in total, with 13 of each suit) provides you with the following (unrealistic) proposition: You have $n$ chances (i.e. $n$ draws with replacement and reshuffling, where $n \to\ \infty$) to randomly select any card from the the suit of spades. If you're successful, you acquire the man's life savings; however, if you're unsuccessful, he acquires yours.* *(In order to eliminate complex scenarios, assume both life savings are equal and that you are risk neutral (i.e. indifferent to risk).) The first interpretation above suggests that the probability of selecting any card from the suit of spades across $n$ draws remains at $p = \frac{13}{52} = \frac{1}{4}$; this implies that the chance if this event occurring is significantly less than the chance of this event **not** occurring ($\frac{1}{4} < \frac{3}{4}$); hence you should not take this bet. The second interpretation suggets that because $n \to \infty$, we should expect the event to occur, on average, once every $4$ times; you have $n$ chances (where $n \to \infty \implies n >> 4$). We should take this bet as the event will (almost?) definitely happen many times across $n$ draws. It should also be noted that the first interpretation suggests that the event may never occur (in fact, it is more likely that it won't occur!), whilst the second interpretation suggests that as $n \to \infty$, the event will occur infinitely many times! I suppose my first question should be: In the scenario above (replacement, independence, etc.), is my reasoning that the probability $p$ of event $E$ occurring across $n$ stagings constant, and, hence, equivalent the probability of the event occurring across one staging of the experiment correct? Secondly, assuming the answer to the previous question is affirmative, is the perceived conflict presented in the above scenario founded on sound principles? If so, how should we progress in the presence of this conflict?


title: How is the complexity of the relationship between pattern and category quantified?

text: Lets say you have a class $C$ of binary strings $s$, and each string corresponds to one of $n$ categories. How do you quantify the complexity of the relationship between the sequence of bits in $s$ and the categorization of $s$. For example, Imagine your strings are 6 bits long, $n=2$, and the first bit is the only one that has any bearing on categorization. That would be pretty simple. On the other hand, you could have cases where the relationship is described by a formula, eg: $\operatorname{Category} =x_5$XOR$x_4$XOR$x_3$XOR$x_2$XOR$x_1$XOR$x_0$ or you could have a case where arbitrary rules indicate the category, eg: "any string of the form -101-- falls into the category 6, unless $x_0$ is true" In any case, how do you quantify the complexity of the correspondence between pattern and category.


title: Why is a deterministic trend process not stationary?

text: I am confused on why a simple trend process is not stationary. Consider the following process: $Y_t = a + bt + \epsilon_t$ The variance is clearly constant. However, the mean $bt$ is dependent on $t$. When shifted in time, the mean only depends on the time interval and is independent of history. For example, $Y_{0,t}$ and $Y_{t,2t}$ would have the same mean and variance. So why is this process not stationary? Secondly, if we now consider the following process: $Y_t = a + \sqrt{t}\epsilon_t$ Assume $\epsilon_t$ is standard normal. In this case, the mean is constant, however, the variance is dependent on $t$. However, in this case the variance is proportional to the time interval, which means $Y_{0,t}$ and $Y_{t,2t}$ would have the same mean and variance. So why is this process not stationary? If you could explain it intuitively rather than definition/proof that would be helpful. My understanding of a stationary process is that the first two moments of the process (mean and variance) remain the same when shifted in time or space. I don't think I have the right understanding of stationary processes.


title: Using R's formula style lm() with SPSS

text: Is there a way to type in a regression formula in SPSS in the same manner as R, with weights? For example, in R, I would write something like:               lm(y ~ B1 + B2 + B1*B2, data=df, weights=x)       How would I go about recreating this in SPSS/PASW?


title: How to increase longer term reproducibility of research (particularly using R and Sweave)

text: **Context:** In response to an earlier question about reproducible research Jake wrote > One problem we discovered when creating our JASA archive was that versions > and defaults of CRAN packages changed. So, in that archive, we also include > the versions of the packages that we used. The vignette based system will > probably break as folks change their packages (not sure how to include extra > packages within the package that is the Compendium). > > Finally, I wonder about what to do when R itself changes. Are there ways to > produce, say, a virtual machine that reproduces the entire computational > environment used for a paper such that the virtual machine is not enormous? **Question:**   * What are good strategies for ensuring that reproducible data analysis is reproducible in the future (say, five, ten, or twenty years after publication)?   * Specifically, what are good strategies for maximising ongoing reproducibility when using Sweave and R? This seems to be related to the issue of ensuring that a reproducible data analysis project will run on someone else's machine with slightly different defaults, packages, etc.


title: Any ideas about how to analyze survival data with pseudo-replication (dependent data)?

text: I am a student, working with a team on a large-scale ecological experiment. We want to analyze survival data which has been derived from an experimental design with some pseudo-replication. This pseudo-replication was not discovered, unfortunately, until the middle of the experiment, at which point the design could not be altered. The experimental design involves comparing the survival data of several treatment groups, each comprised of 10 replicates (aquaria) with 10 individuals in each tank. We have measured mortality as a response variable, in response to different environmental stressors. The trouble is that we cannot say that each of the deaths that occurred, have occurred independently since each tank has many individuals. We would like to acknowledge this problem and address it in our analysis. All of the survival analysis tools we are aware of, assume independence between replicates. We are considering using Kaplan-Meier curves, Cox proportional hazard, or even a glm with Gamma error distributions. Any ideas about how we can properly address this problem in order to detect a difference in survivorship between the treatment groups?


title: Median of Medians calculation

text: I do a bunch of Real Estate Reporting and the Median price is often reported. Particularly by the NAR (National Association Of Realtors) ,,, How can this be,, As best I can figure they only get the Medians from each area,, so how should they calculate a national median, a median of medians a a simple average of medians or weighted average of medians or....??? How valid can it be,, I know that the NAR is not getting the total transaction table so can a reasonable representation of Median (nationally be done? Particularly since regional density and prices and market variances are so large now?


title: How should I analyse cyclic load times?

text: I'm in the process of analysing different PHP frameworks and part of this involves testing the load times. For this, I've programmed various usage scenario's. I'm however quite unsure on how to proceed from this. I'm using Apache Benchmark to simulate the number of requests and the number of concurrent requests. When I put the progression of the load times over time, I get the following chart: ![enter image description here](http://i.stack.imgur.com/WSuE6.png) As you can see, it's cyclic (I'm assuming the server clears a cache of sorts at one point, to drop the load times). When I created a histogram, I got this: ![enter image description here](http://i.stack.imgur.com/gKXZh.png) From looking around here, that seems to resemble a Weibull-distribution. I'm however totally unfamiliar with that and looking into it as we speak. I'm basically looking for a way that let's me compare different sets of similar data to find out which is the fastest. Ideally this would be a single value, but I'm not sure if that's possible. Since I haven't got that much experience with statistical analysis, I was hoping that one of the bright minds here might have the answer (or at least guide me in the right direction). If I've forgotten any information that might be useful, just let me know and I'll try to update as best as I can.


title: How to simulate only stationary AR(1) with φ = 0.9?

text: I am interesting in simulating AR(1) processes with 0.9 parameter and n = 10. The itterations should be 10000. When I was trying to run the program it gave me an error in the estimation procedure. "Error in arima(newyt, order = c(1, 0, 0)) : non-stationary AR part from CSS " My code is this:               set.seed(1123)     rep<-10000     phix<-0.90     t<-c(NA,rep)     tstar<-c(NA,rep)     tdstar<-c(NA,rep)     r<-c(NA,rep)     cut.est<-c(NA,rep)     cut2.est<-c(NA,rep)     for(i in 1:rep)     {     ts.sim3 <- arima.sim(n = 63, list(ar=c(phix)))     new<-ts(ts.sim3[c(0,14:63)])     cut.ts<-arima(new,order = c(1,0,0))     yt<-arima.sim(n = 63, list(ar=c(phix)))     newyt<-ts(yt[c(0,14:63)])     cut2.ts<-arima(newyt,order = c(1,0,0))     r[i]<-cor(new,newyt)     tstar[i]<-r[i]/(sqrt((((1+(phix*phix))/((1-(phix*phix)))/(50)))))     t[i]<-r[i]/(sqrt((1-r[i]^2)/(50-2)))     cut.est[i]<-coef(cut.ts)["ar1"]     cut2.est[i]<-coef(cut2.ts)["ar1"]     tdstar[i]<-r[i]/(sqrt((((1+(cut.est[i]*cut2.est[i]))/     ((1-cut.est[i]*cut2.est[i])))/(50)))))     }     values.t.test<-matrix(t)     mean(cut.est)     mean(cut2.est)     values.t.test     sum(abs(t) > 1.96)     sum(abs(tstar) > 1.96)     sum(abs(tdstar) > 1.96)      Do you have any ideas?


title: Calculating sum of squares between groups

text: I have carried out this ANOVA:               summary(aov(drat ~ cyl, mtcars))                      Df Sum Sq Mean Sq F value   Pr(>F)         cyl          1  4.342   4.342   28.81 8.24e-06 ***     Residuals   30  4.521   0.151                          ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      According to the output, the sum of squares between groups is 4.342. I am trying to calculate the sum of squares using R code First I have taken substracted the means of each group from the values and squared the result:               cyl_4_minus_mean <- (mtcars$drat[mtcars$cyl == 4] - mean(mtcars$drat[mtcars$cyl ==4]))^2     cyl_6_minus_mean <- (mtcars$drat[mtcars$cyl == 6] - mean(mtcars$drat[mtcars$cyl ==6]))^2     cyl_8_minus_mean <- (mtcars$drat[mtcars$cyl == 8] - mean(mtcars$drat[mtcars$cyl ==8]))^2      I then tried to multiply the results by the numbers of values within each group and the sum the result:               sum((11*cyl_4_minus_mean), (7*cyl_6_minus_mean), (14*cyl_8_minus_mean))      This gives 49.4459, which is different to 4.342 given by `summary(aov(drat ~ cyl, mtcars))`. Where am I going wrong here?


title: Trouble minimizing perplexity in LDA

text: I am running LDA from Mark Steyver's MATLAB Topic Modelling toolkit on a few Apache Java open source projects. I have taken care of stop word removal (for e.g. words such Apache, java keywords are marked as stopwords) and tokenization. I find that perplexity on test data always decreases with increasing number of topics. I tried different values of `ALPHA` but no difference. I need to find optimal number of topics and for that perplexity plot should reach a minimum. Please suggest what may be wrong. Definition and details regarding calculation of perplexity of a topic model is explained in this post. **Edit:** I played with hyperparameters alpha and beta and now perplexity seems to reach a minimum. It is not clear to me as to how these hyperparameters affect perplexity. Initially I was plotting results till 200 topics without any success. Now on the same range minimum is reached at around 50-60 topics (which was my intuition) after modifying hyperparameters. Also, as this post notes, you bias optimal number of topics according to specific values of hyperparameters.


title: How is this equation read?

text: I want to understand this paper on brain tumour segmentation. ![enter image description here](http://i.stack.imgur.com/RrCKP.png) How is this equation read? I'm guessing $q_i(t_i)$ represents the likelihood of tumour on voxel i.Is q usually used to represent likelihoods? What does the equal sign with the triangle mean? I understand what each term on the conditional probability represents, but I don't get why it's called _proportional_ to the product of the summation. I'm looking for a reading reference here.


title: Why beta process is useful ?

text: Why Beta process is useful/important? How different is that from Dirichlet process?


title: How to visualize tags?

text: I have a data set which contains a hundred or so user **tags** - that is, binary variables marking whether a user has a certain property or has shown a certain behaviour. (data from a web site.)               user, high_usage, tablet, likes_cats, night_owl     A,    1,          0,      1,          0     B,    0,          1,      0,          1     C,    1,          0,      0,          0      How do I visualize this data? Pairwise or all together. I'd like to create a visualization that allows exploration of this data to find groups of behavior.. Any ideas?


title: "When to use boxplot and when barplot" rules (of thumb?)

text: Both **_box-and-whisker_** plot and **_bar chart_** are appropriate graphics for ANOVA according to The R Book (Crawley, 2013), but **which is more appropriate**? I suppose it depends on situation... can anybody help me?


title: Nonnegativity of model optimism

text: The following sounds like a very basic question in learning theory to me, so I am hoping someone can point out the obvious. Efron's "expected optimism" is the expected difference between the prediction and training errors. Efron (2004) shows that for a wide range of loss functions (the "q-family"), and virtually all modeling approaches, the expected optimism in modeling the mean of a vector is proportional to the sum of covariances of fitted and observed components. It thus seems natural that any sensible modeling approach would have non-negative expected optimism (because one tries at least in part to minimize the training error). I wonder if a more general result exists, though. Suppose that we have an arbitrary optimization criterion (not necessarily a loss in Efron's q-family, maybe not even convex/differentiable/separable; just a function of the data and a tunable parameter). We minimize this criterion on the training data over all possible parameter values in some feasible set. Is the expected minimum obtained this way (where the expectation is taken over all training sets) necessarily an underestimate of the expected criterion value applied to new data, unseen during training (where now the expectation is over all training and all iid testing datasets, as in Efron's definition of optimism)? In other words, is the expected optimism of a sensible modeling approach always negative? Or maybe one has to really require something of the distribution of the data and/or know something about the criterion? Counter examples may also be enlightening.


title: Model broken stick model in R where one line has a constant gradient?

text: I am current using the package (lme4) for a mixed effect model with random effects. My model takes the form:               mod <- lmer(response ~ b1(predict1, bp) + b2(predict2, bp) + (1 | site), data = data)      The functions b1 and b2 come from the helpful advice in: Estimating the break point in a broken stick / piecewise linear model with random effects in R [code and output included]               b1 <- function(x, bp) ifelse(x < bp, bp - x, 0)     b2 <- function(x, bp) ifelse(x < bp, 0, x - bp)      I am wondering, how would I set up my model so that the second part (line) of the broken stick model has a constant gradient? So for example, I would like something like: If X < breakpoint, Y = m1X + C1 If X > breakpoint, Y = C2 Hope I am making sense. Matt.


title: How to select validation data when training a neural network?

text: I am training a neural network with time dependent financial data. In order to avoid overfitting I would like to stop the training at the point where my neural network stops improving on a set of validation data, different from the training data. My question is how to best divide my data in training and validation data? And should validation data be randomly selected or be taken all from the end of the sample? I am trying to predict stock prices based on a window of past prices, and if I randomly select data for validation could this lead to a information leakage?


title: In univariate logistic regression, does the scale of values affect the predicted "risk" at a particular value of the independent variable?

text: I have many datasets representing different populations. Each dataset contains values of an independent variable along with (0/1) representing the occurrence of an undesirable event. I analysed each of the datasets separately by using univariate logistic regression. The goal was to calculate and compare the risk at a particular value in all the datasets. I was then able to put the different populations in order based the risk predicted at the same value V. The problem is that I found that the risk calculated has quite a strong negative correlation with the maximum values in the datasets. In other words, the population X with the lowest risk happens to have the largest maximum in its dataset compared with the other populations. I am now concerned that my conclusion (i.e. the population X has the lowest risk at value V) may only be a result of having the largest independent value in its dataset, and not because of any other factor. Is my concern correct?


title: what does "a distribution over distributions" mean?

text: I am reading a pdf about Dirichlet Process, and it said "A Dirichlet Process is also a distribution over distributions." anyone could explain this in plain English what does it mean by that? thanks a lot!


title: Very Low ICC in HLM

text: I ran a HLM analysis (One Way ANOVA with random effects model). I calculated my ICC and it is very small (00.1%). I understand this means that there is basically no meaningful level-2 variability in the level-1 variable. I have to run an analysis on the data still (it was a hypothesis). Do I still use HLM and just note the low ICC and hence unlikely to find anything of value? I had read before that a significant ICC value means that HLM is warranted. Does this mean that I should use a different analysis (e.g., standard regression)? I am trying to come up with a rationale or explanation if asked a question such as, why did you use HLM with an ICC of 0%?


title: Which metric use for problem of clustering spatial data of wind direction and speed?

text: I have two dimensional spatial (x,y - coordinates of meteo stations) data for small region (so I could neglect the shape of earth globe), for each (x,y) I have one observation of wind direction and speed I need to design metric for clustering problem of such data. It's simple to do when we consider only directions or only the speed of wind but how to combine these quantities ? I try to analyse this problem using R.


title: Shifted intercepts in logistic regression

text: I have a question about the effects of shifting the intercept in a logistic fit on the mean of a particular transformation of the scores. Here is the notation I will be using for the question. The logistic regression model is $$ \begin{align} \mathbb{P}[Y_i =1 \mid \boldsymbol{X}_i] &= \dfrac{1}{1+\exp(-\boldsymbol{X}_i'\boldsymbol{\beta})} \\[1em] &\equiv p(\boldsymbol{X}_i; \boldsymbol{\beta}) \end{align} $$ where $p$ is the score function. I estimate the parameters of the model, and denote them $\widehat{\boldsymbol{\beta}}$, so that the estimated scores of the model are $p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})$. I am interested in the sample mean of all the estimated scores that cross a certain threshold, $\delta$, say. That is, the quantity $$ \dfrac{\sum_{i=1}^n p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\bf{1}_{\left[p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]}}{\sum_{i=1}^n\bf{1}_{\left[p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]}} $$ Now suppose that I shift the estimated linear index for all observations by the same quantity $\alpha$, and denote the transformed scores as $$ \begin{align} q(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}}) &= \dfrac{1}{1+\exp(-(\boldsymbol{X}_i'\widehat{\boldsymbol{\beta}}+\alpha))} \\ &= \dfrac{\exp (\alpha)}{\exp (\alpha)+\exp(-\boldsymbol{X}_i'\widehat{\boldsymbol{\beta}})} \end{align} $$ It is easy to see that if $$ \begin{alignat}{2} & \alpha &\geq& 0\\ \implies\,& \exp(\alpha) &\geq& 1\\ \implies\,& q(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}}) &\geq& p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}}) \\ \implies\,& \sum_{i=1}^n q(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\bf{1}_{\left[q(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]} &\geq& \sum_{i=1}^n p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\bf{1}_{\left[p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]}\,\text{and also,}\\ \implies & \sum_{i=1}^n\bf{1}_{\left[q(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]} &\geq& \sum_{i=1}^n\bf{1}_{\left[p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]} \end{alignat} $$ However, since the denominator and the numerator both become larger for the $q$ function, it is hard to tell what will happen to the ratio, $$ \dfrac{\sum_{i=1}^n p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\bf{1}_{\left[p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]}}{\sum_{i=1}^n\bf{1}_{\left[p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]}} $$ that is $$ \dfrac{\sum_{i=1}^n q(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\bf{1}_{\left[q(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]}}{\sum_{i=1}^n\bf{1}_{\left[q(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]}}\overset{?}{\lesseqqgtr} \dfrac{\sum_{i=1}^n p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\bf{1}_{\left[p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]}}{\sum_{i=1}^n\bf{1}_{\left[p(\boldsymbol{X}_i; \widehat{\boldsymbol{\beta}})\geq \delta\right]}}$$ I am pretty sure there is an example like this in _Theory of Point Estimation_ by Lehmann and Casella, for the univariate (non-regression case), but I have been unable to find it. * * * # Simulations I decided to simulate the above scenario to see if I could get an intuition for this problem. The data generating process that I have used for generating the data is slightly non-standard, since I wanted to make sure that my simulations were robust to distributional assumptions. It might be useful to try different parameters configurations in choosing the fixed parameters. Here is the `R` code for the simulations (I would appreciate it if someone could confirm that this code is indeed doing what I think it is).               # purpose: simulate logistic regression     # date created: 27th October 2012     # comments:     # author:     # dependencies: faraway          library(faraway)          # simulation parameters     vBeta <- rbind(0.1, 0.2, 0.3, 0.4)  # vector of coefficients     sDelta <- 0.16  # threshold for the scores     sAlpha <- 0.7  # shift in linear index          # matrix to hold the results     mResults <- matrix(0, nrow=1000, ncol=2)          for(i in seq(from=1, to=1000, by=1)) {       # simulate the data       mX <- matrix(rnorm(400, 4, 1), 100, 4)       vY <- (0.4 + mX%*%vBeta + rt(n=100, df=7)>=5)       data <- as.data.frame(cbind(vY,mX))            # logistic regression       resLogitFit <- glm(V1~V2+V3+V4+V5, binomial(link = "logit"), data=data)            # save the mean of the scores that are above a certain value       mResults[i,1] <- mean(resLogitFit$fitted.values[resLogitFit$fitted.values > sDelta])            # effect the logit transformation on the shifted linear index       data$shiftedScores <- ilogit(resLogitFit$linear.predictors+sAlpha)       mResults[i, 2] <- mean(data$shiftedScores[data$shiftedScores > sDelta])     }          # plot the two means across simulations     matplot(y=mResults, type="l",ylim=c(0.2, 1), xlab="Simulation #",              ylab="Estimated statistic")     legend(1, 1, c("(mean) truncated scores", "(mean) shifted truncated scores"),            cex=0.8, col=c("black","red"), lty=1:2)          # count the number of times the shifted scores exceed the older means     sum(mResults[,1] > mResults[,2])      This produces the following graph, which seems to indicate that the shifted truncated scores have a higher mean -- although this depends on the shift parameter. ![enter image description here](http://i.stack.imgur.com/Gkom2.png) I would very much like a theoretical explanation to support these simulation results.


title: Joint probability and Gaussian copula

text: I have $\Pr(A)=29\%$ and $\Pr(B)=10\%$, where $A$ and $B$ are two events which are not independent. In fact, a correlation measure suggests they're correlated by $\rho=0.8$. I would like to calculate $\Pr(A \cap B)$, that is the probability both of them occur. Here is what I've done in R:               prA <- .29     prB <- .1     Sigma <- matrix(.8, 2, 2)     diag(Sigma) <- 1     qA <- qnorm(prA, 0, 1)     qB <- qnorm(prB, 0, 1)     library(mnormt)     pmnorm(x = c(qA, qB), mean = rep(0, 2), varcov = Sigma)      Is the result of such a code $\Pr(A \cap B)$ I'm looking for?


title: Weighted Distributions in Predictive Model

text: I'm not a statistician but do work with large datasets and have a problem I'd like to use a predictive model for. I have two datasets that I'd like to use together to build predictions. The first set is of a list of categories: A, B, B, A, C, D, A, B, D , etc. I have turned this into a table of histograms. Essentially, I have created categorical probabilities for each of the labels above. **Let's call this distribution X.** The second set is a table with counts by source and category, here is a very simple example, though there are many more for each source:               source   cat   count     one      A     19     one      B      0     one      C      10     one      D      1     two      A      0     two      B     20     two      C      1     two      D      0     three    A     100     three    B      30     three    C      57     three    D      3      I've also created a distribution for each source, call this **distribution Y**. What I'd like to do is given cat A or cat B as a possible next entry, calculate which is more likely and with what margin of error. To do this I will use both the X and Y with a weight **for each source** : Predicted_prob(cat A) = w * X(cat A) + (1 - w) * Y(cat A) Here are some questions about this approach:   1. Is this an ok general approach? I don't need anything super fancy, just something that can work with the data I have and give a way of measuring the confidence of the prediction.   2. I am currently using a uniform distribution for w where w = (1 / n) and n is the number of samples of each source. Is there a better distribution that may work to make sure if n is too small, the effects are not felt as much?  I am working out the research in R but will ultimately need to implement with Python (can use RPy). I'm open to any theoretical or implementation input, as long as it's attainable!


title: Jonckheere-Terpstra interpretation

text: I am running the Jonckheere-Terpstra in place of Kruskal-Wallis test, as my factor is in ordinal scale (i.e. groups can be ordered). The Asymptotic significance (2-tailed) is 0.000, so it seems there is a trend in the response variable, according to the trend of the factor. However, if I look at the medians, this trend is not clear; in fact I have 4 levels in my factor and, in ascending order, their medians are: 0.1387, 0.2814, 0.5882, 0.3492. So, I ran a Mann Withney test with Bonferroni correction to check if the difference between last two pairs is significant and it is NOT. Two questions: 1) Am I doing it right? 2) I ran the J-T test for another response variable, which is ordered conversely with respect to the factor (i.e. if I increase the factor, the response variable decreases). If the J-T statistic is NEGATIVE and the p-value is 0.000, can I conclude that the response variable is DECREASING (instead of increasing) according to a factor increase? Thanks.


title: How can I learn about fitting nonlinear and nonadditive models?

text: The classical (OLS) regression model is: $$ y_i = \beta_1 x_{i1} + \cdots + \beta_p x_{ip} + \varepsilon_i $$ For this linear model we have many instruments to estimate model's fit. But how can I estimate the fit of model like this: $$ y_i = \beta_1 x_{i1} * \cdots * \beta_p x_{ip} * \varepsilon_i $$ or other model with division / log, etc? I want to find information about methods for this task (keywords, references, articles, etc)


title: How to do external validation of logistic regression models and perform model benchmarking

text: Quality assessment in trauma has for > 25 years been done with the US derived logistic regression model, the TRISS model. DV: survival/death and IVs: physiologic derangement (continuous), anatomic injury (continuous), age (dichotomous, >= 55, <55). Probability of survival (Ps) are calculated for each patient: $$ Ps = \frac{\exp\big(b_0 + b_1{\rm (physiology)} + b_2{\rm (anatomy)} + b_3{\rm (age)}\big)}{1+\exp\big(b_0 + b_1{\rm (physiology)} + b_2{\rm (anatomy)} + b_3{\rm (age)}\big)} $$ Used worldwide with the US derived regression coefficients, updated in 2005 and in 2009, it has had huge impact in trauma research. Foreign institutions has been able to benchmark their performance towards the US standard by use of the W-statistic (with 95% CI), expressing excess (less) survivors per 100 patients in own institution vs. the US standard. (W = (actual number of survivors – predicted number of survivors)/(number of patients/100). Some countries have used the same DV and IVs, however derived their own regression coefficients from their own population. Others have derived their own logistic regression model. We have recently derived our own logistic regression model with the same DV, but different IVs. We have also implemented the TRISS model in our registry with the national US IV, but with IV coefficients derived from our trauma population. Questions:   1. I want to perform an external validation of the US model with national US IV regression coefficients. I want to compare the predictive ability/performance vs. our own model. How should we perform such a comparison? Compare ROC curves, deviance, or AIC? Is it a meaningful comparison?    2. I want to perform an external validation of the US model, however with IV regression coefficients derived from our trauma population. I want to compare the predictive ability/performance vs. our own model. How should we perform such a comparison? Compare ROC curves, deviance, or AIC? Is this a more relevant comparison? How can we tell whether differences are significant or not?    3. What about using Net Reclassification Improvement?    4. I want to decide which model fits our data best. How?    5. Can other institutions use our model for benchmarking?


title: Time series image data classification / video image classification

text: I am working on classifying video frames into two classes, positive and negative. e.g. if a particular pattern appears in a frame that frame will be classified into positive, otherwise negative. But in some cases it is not clear whether that pattern appears in a particular frame (let's say i) or not. But observing the adjacent frames of i, I can easily classify whether i is positive or negative. I am not sure how to use this temporal information to classify a particular frame. May be I can use Kalman filters (KF) or Markov random fields (MRF), etc. But I am unfamiliar with KF or MRF. So could you please give me some suggestions and some guidance in the recent literature to accomplish this?


title: Dickey Fuller test

text: I am studying the Dickey Fuller test. The book of reference is Introductory Econometrics for Finance by C. Brooks. I firstly consider the zero mean Dickey Fuller test that uses the "random walk" type of regression: $\Delta Y_t = \delta Y_{t-1} + u_t$ The null hypothesis is that $\delta$ is zero and the alternative is that $\delta$ is not zero so I run the regression and compare the t-value with the critical values tabulated by Dickey and Fuller. Though the results are valid only if $u_t$ is white noise. According to the book previously mentioned, if the dependent variable of the regression (the first difference of the series) is autocorrelated then the $u_t$ would also be autocorrelated. Why is that, how this can be proved? Also if the dependent variable is not autocorrrelated can we infer that the $u_t$ of the regression would also be uncorrelated? In general can we draw coclusions about whether the error term is autocorrelated or not based on whether the dependent variable is autocorrelated? If yes, can the same inferenece be made in the case of the two other forms of the Dickey Fuller test regressions (1) the equation above with a constant term (drift) or single mean, 2) The equation above with drift and time trend).


title: Using a nominal independent variable in scikit-learn

text: I am trying to understand how scikit-learn uses input data to train and build a classifier. Please note that i am both a Python and scikit-learn newbie. From what i understand so far, scikit-learn does not handle categorical features as with other Software (IMHO a serious drawback if this actually holds) such as R, WEKA, KNIME, RapidMiner, IBM SPSS etc. The data contains 5 numerical features and the target nominal variable (values TRUE/FALSE) I managed to go through the examples on scikit learn and came up with the following code :               from numpy import genfromtxt,zeros               data=genfromtxt('scores.csv',delimiter=',',usecols=(0,1,2,3,4),skip_header=1)     target=genfromtxt('scores.csv',delimiter=',',usecols=(5),dtype=str,skip_header=1)          print data.shape          print set(target)               t=zeros(len(target))     t[target=='TRUE']=0     t[target=='FALSE']=1               from sklearn import tree               classifier = tree.DecisionTreeClassifier()     classifier.fit(data,t)                from sklearn import cross_validation     train, test, t_train, t_test =cross_validation.train_test_split(data,t,test_size=0.4, random_state=0)          classifier.fit(train,t_train) # train     print classifier.score(test,t_test) # test          from sklearn.metrics import confusion_matrix     print confusion_matrix(classifier.predict(test),t_test)          from sklearn.metrics import classification_report     print classification_report(classifier.predict(test), t_test, target_names=['FALSE', 'TRUE'])               from sklearn.externals.six import StringIO       with open("iris.dot", 'w') as f:     f = tree.export_graphviz(classifier, out_file=f)      The tree output (and more specifically the leaves)produced by the algorithm does not contain the original class labels but output such as : value = [ 127. 0.] Is it possible to have the actual class labels shown rather than numbers?


title: Data sets and problems for learning hypothesis testing

text: I'd like to get some real practice with hypothesis testing that is above what's in textbooks (plug and chug, usually). I was hoping someone could suggest some good data sets and problems to work on. Problems that I could show others that I've worked on. Along those lines, I'd very much appreciate any top notch references on the topic.


title: How to correctly describe measurement variation (over multiple repetitions)

text: Firstly, please forgive my lack of statistics knowledge, but I am hoping someone can clear up my misunderstandings. I am taking a sample of a physical quantity, (for example: temperature). To take this sample I average the measurements over a time duration (for example: 1 sec). This gives me a sample mean and a sample deviation. I then repeat this measurement procedure under the same conditions. I then would have a sample mean and sample standard deviation for each repetition. So my question is, while each sample has variation how do I "combine" these to describe the variation over all of the repetitions? This seems rather fundamental, so please link any relevent references if you know of any.


title: How to find 1-year population average with intermittent series data

text: I have performance review data and scores (ranging from 1 to 4) for employees of a company. I need to show the company average over the past year. However, the employees were only ever reviewed for a few weeks in a row at a time. For example:   Roger was reviewed from Jan 1 through March 1   Steve was reviewed from August 15 through September 15   ...etc... Right now, I'm taking all of the reviews over the year (Jan 1 through Dec 31) and for each day I'm computing the average for that day. I then take all of those days and plot the average score for that day on a line graph with time on the X axis and score on the Y axis. The problems I'm having are:   * there are gaps of time in these reviews (no one was reviewed in February, 20 people were reviewed in March)   * the higher performing employees may have been reviewed at one part of the year, and the lower performing employees at another   * when I plot this data the average produces a zigzag looking line chart. for example, day 1 the average is 3.4, but day 2 is 1.3. over time this looks like a lot of spikes. So, where I'm at now is thinking of taking each employee as a "series" (where some series (employees) have data ranging for a few weeks, and some have data ranging a few months) and normalizing them together and computing the average of that normalization. I'd greatly appreciate any help on this and the best way to present this data!


title: Conceptual links and differences between generalized and latent models

text: Please bear with me, as I'm trying to better understand these aspects. My understanding is that _(finite) mixture models (MM)_ are characterized by a presence of a number of sub-populations in a population or sample. It can be argued that this _heterogeneity_ might represent effects of some known or unknown _factors_ , internal or external. This sounds related to the concept of _generalized additive models (GAM)_ , which are concerned with some unknown _predictor variables (factors)_. On the other hand, we have such terms as _latent Gaussian models (LGM)_ as well as _latent variable models (LVM)_ , which, as I understand, is just another term for _structural equation models (SEM)_. Assuming that this _logic_ is conceptually (in a big picture way) correct, my **question** is two-fold: 1) what are major differences and relationship _between MM and GAM_ concepts; 2) what are major differences and relationship _between MM/GAM and LGM/LVM_? In other words, what are concepts, theories and approaches **connecting** _SEM_ and other types of models mentioned by me above? Or, turning this question on its practical side: to what extent and about which aspects do researchers need to care about all these models, when their main approach in a particular study is SEM?


title: Adaptive stopping rule for spatially autocorrelated Poisson data

text: I have a problem where I am given an initial prior or proportion for the number of occurrences per some unit with known standard deviation (an example is 3 per mile). I wish to test if this proportion is valid given a live exhaustive sampling of the area. So I would plan to use all the sampled information. I understand that something like a z test for proportions works when you are working with a given sample size (as this is used in calculating a test deviation). One way I have though of solving this is to simply break down the live data into a a series of small 1 meter size samples. In addition the independence assumption may be broken because the sampling is not truly random given the live data is likely going to be moving from area to area. That is, at any given point there will be large areas of the domain that have never been sampled from until the end. So the question becomes: **is there a way to use a continuous sample to perform a hypothesis test, and is there any way to account for the biases in the sampling process?**


title: How do I generate numbers according to a Soliton distribution?

text: The Soliton distribution is a discrete probability distribution over a set $\{1,\dots, N\}$ with the probability mass function $$ p(1)=\frac{1}{N},\qquad p(k)=\frac{1}{k(k-1)}\quad\text{for }k\in\{2,\dots, N\} $$ I'd like to use it as part of an implementation of an LT code, ideally in Python where a uniform random number generator is available.


title: Simulating a price series using GARCH

text: I need to simulate a price series based on a GARCH(1,1) specification for the returns (price changes). I currently have this:                   d_price = diff(price) #the price changes         garch.model = garch(d_price) #a GARCH(1,1) model         specs = garchSpec(model = list(omega = garch.model$coef[1],              alpha = garch.model$coef[2], beta = garch.model$coef[3]))              sim = garchSim(spec = specs, n = t) #simulated GARCH returns         mu = mean(price) #the starting value is the mean of the series         sim.series = rep(mu, t)               #create the simulated price series         for(i in 2:t){            sim.series[i] = sim.series[i-1] + sim[i-1]         }      My problem is that this series often becomes negative. I can't have negative prices in my model. I thought about using percentage price changes instead, but I experience a similar problem - the price series sometimes decreases to arbitrarily small values that are completely unrealistic. **So, my question is this:** Is there a best option for simulating a series so that the price never becomes negative or unrealistically small? That is, without putting artificial barriers on the price process (P must be greater than 20, for example). Thanks!!!


title: Supervised clustering or classification?

text: The second question is that I found in a discussion somewhere on the web talking about "supervised clustering", as far as I know, clustering is unsupervised, so what is exactly the meaning behind "supervised clustering" ? What is the difference with respect to "classification" ? There are many links talking about that: http://www.cs.uh.edu/docs/cosc/technical-reports/2005/05_10.pdf http://books.nips.cc/papers/files/nips23/NIPS2010_0427.pdf http://engr.case.edu/ray_soumya/mlrg/supervised_clustering_finley_joachims_icml05.pdf http://www.public.asu.edu/~kvanlehn/Stringent/PDF/05CICL_UP_DB_PWJ_KVL.pdf http://www.machinelearning.org/proceedings/icml2007/papers/366.pdf http://www.cs.cornell.edu/~tomf/publications/supervised_kmeans-08.pdf http://jmlr.csail.mit.edu/papers/volume6/daume05a/daume05a.pdf etc ...


title: Goodness of fit for stock prices

text: I'm in my first days of learning `R` and hit a roadblock with a small use case from finance. * * * (EDIT) Basically I want to know how to test sample fitness against any distribution. Say, find the degree of freedom of a t-distribution which would have the best fit with the given list of asset returns. I thought the way to go is to generate a random dataset with given properties, and `qqplot` it against the actual data, along with an `abline` to see if it's linear or not. In the question I'm using normal/lognormal just as an example which I thought would be easier to discuss. * * * So, there's a list of stock prices (presumably a lognormal distribution), from which I calculate a list of returns (presumably a normal distribution). Now I want to check my assumptions about the distributions, so I plot both:               par(mfcol=c(2,3))          #plot densities     plot(density(sbux.prices))     plot(density(sbux.returns))          #plot q-q probabilities for normal distribution     sbux.prices.norm = rnorm(n=1000, mean=mean(sbux.prices), sd=sd(sbux.prices))     qqplot(sbux.prices.norm, sbux.prices)     abline(0,1)      sbux.returns.norm = rnorm(n=1000, mean=mean(sbux.returns), sd=sd(sbux.returns))     qqplot(sbux.returns.norm, sbux.returns)     abline(0,1)          #plot q-q probabilities for lognormal distribution     sbux.prices.lnorm = rlnorm(n=1000, mean=mean(sbux.prices), sd=sd(sbux.prices))     qqplot(sbux.prices.lnorm, sbux.prices)     abline(0,1)      sbux.returns.lnorm = rlnorm(n=1000, mean=mean(sbux.returns), sd=sd(sbux.returns))     qqplot(sbux.returns.lnorm, sbux.returns)     abline(-1,1)      ![enter image description here](http://i.stack.imgur.com/yDix8.png) Two questions here.   1. By just looking at the plots above, I can tell that prices follow lognormal distribution - because the corresponding QQ plot has a much better fit with lognormal than with normal. But returns seem to fit well with both distributions - is that correct or I'm doing something wrong?    2. Obviously I'd better rely on a mathematical estimation of fitness rather than on a human checking the charts. I think, I can use Chi-squared or Kolmogorov-Smirnov test, but cannot understand how exactly to do that. E.g., `ks.test(sbux.returns, sbux.returns.norm)` gives me `p-value = 0.007781` so I'm definitely missing something.


title: Can causality be inferred in a study with an experience followed by two sets of measures

text: I came across this study as part of a mock exam paper and was confused to say the least. ### Context: The study investigates cognitive and behavioural factor related to the experience of anxiety in MRI scanners. Participants completed the following questionnaires 5 mins after the scan:   * a measure of the frequency of anxiety experienced in the scanner   * a measure of the frequency of physical symptoms of panic   * a measure of the frequency of coping strategy employment   * a measure of the frequency of claustrophobic-related thoughts Participants were also asked to fill in and submit general measures at home after the scan:   * level of claustrophobic fears regarding restriction and suffocation   * depression   * anxiety   * health anxiety The return rate was 75%, i.e. only 97 of the original 130 participants completed the second batch of questionnaires. Overall, the level of anxiety in the scanner was RELATED TO:   * level of claustrophobic fears regarding restriction and suffocation   * frequency of coping strategy to manage scanner anxiety   * no. symptoms of panic during the scan   * no. claustrophobia-related thoughts about the scan. No stats are given. ### Questions:   * What type of design is this?   * Is it possible to understand the causes of anxiety experienced in MRI using this type of design?  ### Initial Thoughts: In my mind, causality can not be inferred because there is no comparison with a control group, and correlation would probably have been used rather than AVOVA. Is this the answer to the questions above or have I missed the point?


title: Difference between OLS and FE

text: I have a more general question. Could somebody please explain what is the general difference between OLS and FE (Fixed Effects) in a very simple way? In terms of use in panel data and in general. Thank you!


title: Regression analysis (OLS) in matched case-control study

text: I have been asked to assist in a study based on samples of records from registers of dental care patients. My problem is that I would like to be sure that my proposed approach for the analysis is sufficiently common, but due to my own limited experience in the field, I do need advice. One aspect is the final analysis:   * For each case there is one control matched on sex and age.    * It is not a question of calculating odds ratios for any disease or similar as is often done in case-control studies, but rather to estimate differences in the mean number of teeth. So far, I have thought of making the analysis in one of the following ways:   * Calculate the differences between the number of teeth for each case and their controls, i e more as paired-sample comparison, but put into a ordinary regression context making it possible to adjust for sex and age;   * Use the cases number of teeth as the dependent variable and the corresponding control's number as an indep variable. Is this is reasonable or not, are there better alternatives? Any suggestions and comments are welcome.


title: Is there an anova procedure that doesn't assume equal variance?

text: When conducting a t-test we can always use a Welch approximation of the df, and let go of the assumption of equal variances. Is there something similar for (even one way) anova? (any R implementation of that?)


title: Accommodating entrenched views of p-values

text: Sometimes in reports I include a disclaimer about the p-values and other inferential statistics I've provided. I say that since the sample wasn't random, then such statistics would not strictly apply. My specific wording is usually given in a footnote: > "While, strictly speaking, inferential statistics are only applicable in the > context of random sampling, we follow convention in reporting significance > levels and/or confidence intervals as convenient yardsticks even for > nonrandom samples. See Michael Oakes's _Statistical inference: A commentary > for the social and behavioural sciences_ (NY: Wiley, 1986). On a couple of occasions--once for a peer-reviewed paper, once or twice in a non-academic setting--the editor or reviewer objected to this disclaimer, calling it confusing, and felt that the inferential findings should simply stand as written (and be given the mantle of authority). Has anyone else encountered this problem and found a good solution? On the one hand, people's understanding of p-values is generally dismal, even in the context of random sampling, so perhaps it doesn't matter much what we say. On the other, to contribute further to misunderstandings seems to make one part of the problem. I should add that I frequently deal with survey studies, where random assignment does not apply and where Monte Carlo simulations would often fail to address the issue of representativeness.


title: Testing "trends over time" of dummy variables

text: Let us say this is an output of a model I ran in Stata, where `int_ret`is a continuous variable and `time1`-`time17` are time dummy variables. Reference is `time0`. I notice that the coefficients are significant at alpha=0.05 up til time 5 after which it stops being significant. I also notice that the coefficients for time decrease in magnitude from time1 to time5. I am wondering if it is allowable for me to test if there was a "decreasing trend" between time1 and time 5 using the `test` command in Stata. Example: `test time1=time2=time3=time4=time5=0` or                 test time6=time7=time8=time9=time10=time11=time2=time13=0      I know that you have to enter all the dummy variables (except reference) in a model, but can you test the significance **just a subset** of the dummy variables (meaning, NOT all dummy variables) to determine "trend?" Thanks                        |             Linearized      int_ret |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]        time1 |   3.803034   1.569459     3.24   0.001     1.690201    8.557009        time2 |   2.983938   1.145295     2.85   0.005     1.403566    6.343762        time3 |   2.963121   1.412266     2.28   0.023      1.16145    7.559592        time4 |   2.709894   1.137635     2.37   0.018     1.187634    6.183323        time5 |   2.319515    1.00564     1.94   0.053     .9894511    5.437508        time6 |   3.240125   2.107322     1.81   0.071     .9026434    11.63074        time7 |   1.474064   .6102438     0.94   0.349     .6534528    3.325207        time8 |   1.348993   .6327848     0.64   0.524     .5366498    3.391004        time9 |   1.211208   .4933164     0.47   0.638     .5440385    2.696546       time10 |   .7559254   .4048442    -0.52   0.602     .2638907    2.165379       time11 |   2.224805   1.034081     1.72   0.086     .8925419    5.545686       time12 |   1.676299   .6908775     1.25   0.211     .7457992    3.767742       time13 |   1.233598   .6585384     0.39   0.694     .4321069    3.521729


title: Advice needed on optimal sensor selection and rules inference

text: I observe a certain system with a rectangular array of sensors, that can be either triggered or not. I know the system behavior can be well described by the sequence of triggering events and by the time difference between them. The system is not deterministic so it cannot be fully captured by a finite state machine or something similar. Based on a certain family of sequences, the system behavior can be labeled in many discrete ways such as "good", "bad”, “neutral", etc. I also know that probably only a few of the sensors are conveying the really important information. Given the above, I would like to start from a set of sample sequences with known labels in order to:   * Learn what is the subset of sensors that best explain the system behavior in order to start with the complete array, and then find a smaller group of them that still can describe effectively the system.    * Define the group of rules that best capture a certain label.  A rule would be something like: > if S1, S3, S2 with S3 within 1 second from S1 then the system is yellow The rules should be something like firewall rules, in which one can later add some rules that the algorithm learned later after additional training. **What could be a suitable machine learning algorithm or framework to capture the above?** If it is important, the numbers at play are:   * 200-300 sensors that could be more coarsely grained down to 100 if mathematically intractable    * a dozen typical labels    * a training set of 500-1000 learning sequences each made by 200-300 on/off time stamped events


title: How to compare whether models built using two different outcomes are significantly different

text: I would like to build the relationship between the dose given and two outcomes (one acute toxicity and one late toxicity). The model I used was binary logistic regression. For the acute toxicity outcome, the estimated dosage and its 95% confidence interval was 20 (15~35), while for the late toxicity outcome, the estimated dosage and its 95% confidence interval was 75 (50~80). We can see that 75>20 and these confidence intervals are not overlapping. It seems to indicate that the late toxicity is caused by a higher dose than the dose for acute toxicity. My question is: Is there any statistical method to test that the late toxicity is caused by a higher dose than acute toxicity? Or is simply claiming that their CIs do not overlap is strong enough? Thank you


title: KL distance in R

text: I am trying to compute the KL-distance between two columns in a data-frame. Whats the best way to implement this? I have seen FNN and FLEXMIX packages and their examples, but it is implemented for random samples. How can I do it on two data frame columns?                    V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11     1 1000025  5  1  1  1  2  1  3  1   1   2     2 1002945  5  4  4  5  7 10  3  2   1   2     3 1015425  3  1  1  1  2  2  3  1   1   2     4 1016277  6  8  8  1  3  4  3  7   1   2     5 1017023  4  1  1  3  2  1  3  1   1   2     6 1017122  8 10 10  8  7 10  9  7   1   4      My data looks like as above. Now I am trying to computer KL distance between {V2,V3},(V2,V4}... so on. **What did I try till now?**               install.packages("FNN");          library("FNN");          > X <- rexp(BreastCancerData$V2,rate=0.2)         > Y <- rexp(BreastCancerData$V3,rate=0.2)     > KL.dist(X,Y,k=1)     [1] 0.1391855      Is this a right way to do it?


title: How do I estimate a factorial design with empty factors?

text: I'm trying to analyse data from a factorial design, but I don't have all combinations of factor levels, and I have repetition for the combination. As this configuration differs from book's examples I'm wondering how it would impact the results. An example of what I'm saying:               F1 2 levels       F2 3 levels        usually I do               F1.1,F2.1       F1.1,F2.2       F1.1,F2.3       ...       F1.2,F2.3        instead I have               F1.1,F2.1       F1.1,F2.1 (repetition)       F1.1,F2.2 (this combination missing)       F1.1,F2.3       F1.1,F2.3 (repetition)       ...       F1.2,F2.3       F1.2,F2.3 (repetition)


title: Product of Gamma by Beta rv

text: If $X$ has a beta distribution $ \beta(\alpha,b)$, $Y$ has a gamma distribution $\Gamma (K,\theta)$ and $X$ is independent of $Y$. What is the distribution of the product $P=XY$ . Thanks!


title: How do you use simple exponential smoothing in R?

text: I'm beginner in R, Could you please explain how to use ses in forecast package of R forecast? I'd like to choose the number of initial periods and smoothing constant.               d=[3    4   41  10  9   86  56  20  18  36  24  59  82  51  31  29  13  7   26  19  20  103 141 145 24  99  40  51  72  58  94  78  11  15  17  53  44  34  12  15  32  14  15  26  75  110 56  43  19  17  33  26  40  42  18  24  69  18  18  25  86  106 104 35  43  12  4   20  16  8     ]      I have 70 periods, I'd like to use 40 Periods for initial and 30 for out-of sample.               ses(d, h=30, level=c(80,95), fan=FALSE,initial=c("simple"), alpha=.1)      Is it correct?


title: Z-scores and normalization in previously transformed data sets

text: I have two sets of data. The first is multivariate linear regression data consisting of betas and standard errors. It was log(e) transformed before the linear regression model was applied. The second set is also multivariate linear regression data, consisting of betas and standard errors, but it was z-score transformed before the model was applied. I would like to compare these two sets of data, but unfortunately, I do not have access to the raw data and therefore cannot transform them in the same way, from scratch. I would like to apply z-score transform the first set, using the betas, but am wondering:   1. if this makes sense, mathematically,    2. if the z-score formula should change at all to account for the application for the regression, and    3. if there are any other options or considerations I should be reviewing?